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.  On  the  other  hand,  ordinary  tomography  is  not  doing  well  in 
explaining  the  amplitudes,  whereas  diffraction  tomography  explains  about  20%  even 
when  amplitude  data  are  not  used  in  the  inversion,  and  about  40%  when  the  inversion  does 
include  these  data.  It  is  concluded  that  the  modified  method  of  diffraction  tomography 
described  in  this  paragraph  promises  to  become  an  important  tool  for  determining  the 
structure  beneath  seismic  arrays. 

In  paragraph  2.2,  the  problem  of  discriminating  between  earthquakes  and  underground 
explosions  is  formulated  as  an  exercise  in  pattern  recognition  approach  analysis.  The  anal¬ 
ysis  has  been  applied  to  a  learning  set  of  44  nuclear  explosions  (8  test  sites)  and  35  earth¬ 
quakes  in  Eurasia  recorded  at  the  NORESS  array.  The  signal  features  considered  were  the 
normalized  power  in  8  spectral  bands  in  the  0.2-10  Hz  range  of  the  P-wave  (6  sec)  and  the 
P-coda  (30  sec).  Brysically,  it  means  that  we  exploit  potential  differences  in  the  shape  of 
earthquake  and  explosion  spectra,  respectively.  Other  features  included  are  peak  P  and  P- 
coda  amplitude  frequencies  and  relative  P/P-coda  power.  These  19  features  were  extracted 
either  from  conventional  array  beam  traces  or  the  optimum  group  filtered  traces  (GGF- 
removal  of  coherent  low-frequency  noise).  Using  the  feature  selection  algorithm,  based  on 
estimates  of  the  expected  probability  of  mi.sclassiiication  (EPMC),  only  2  to  4  features 
were  needed  for  optimum  discrimination  performance.  The  dominant  features  were  coda 
excitation  and  P-  and  P-coda  power  at  lower  signal  frequencies.  Furthermore,  feature 
parameters  extracted  from  the  OGF  traces  had  a  slightly  better  performance  in  comparison 
to  those  extracted  from  beam  traces.  Finally,  there  were  no  misclassihcations  for  OGF- 
derived  features  when  the  explo.sion  population  was  limited  to  E.  Kazakh,  while  including 
events  from  the  other  test  sites  lead  to  a  decrease  in  discrimination  power. 

We  have  previously  applied  the  ray  tracing  technique  of  Bo.stock  anJ  Kcnnett  ( 1990)  to 
the  characterization  of  Lg  propagation  to  the  NORESS  array.  In  paragraph  2.3,  we  have 
considered  the  application  of  the  same  approach  of  representing  the  major  features  of  the 
Lg  phase  via  constructively  interfering  S  waves  multiply  reflected  within  the  earth’s  ernst 
to  understanding  azimuth  anomai'es  observed  at  the  ARCESS  array.  The  motivation  for 
the  theoretical  modelling  has  been  the  rapid  changes  in  both  size  and  sign  uf  ARCESS  L;; 
azimuth  anomalies  with  geographical  position,  clearly  observed  in  IMS  analysis  of  events 
to  the  south  and  east  of  Rnland  (Gulf  of  Estonia  -  Leningrad  •  USSR  north  of  Leningrac ). 
Theoretical  azimuth  anomalies  calculated  from  this  ray  tracing  scheme,  and  based  on  r 
simplifled  model  of  the  Moho  structure  under  the  Fennoscandian  region,  are  found  *o  i)e 
generally  consistent  with  the  trends  in  the  observed  anomalies  at  ARCESS.  However  the 
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ray  picture  for  Lg  is  probably  more  valuable  for  providing  a  tool  to  investigate  the  sensi¬ 
tivity  of  guided  wave  propagation  to  the  relative  position  of  the  source  and  major  subsur¬ 
face  features  such  as  the  significant  Moho  topgraphy  under  Finland. 

Paragraph  2.4  is  a  contribution  on  the  use  of  Lg  as  a  yield  estimation  tool.  In  recent  years, 
the  Lg  phasve  has  emerged  as  maybe  the  most  promising  tool  to  obtain  precise  yield  esti¬ 
mates  by  seismic  means.  The  pioneering  work  by  Otto  Nuttli  of  Saint  Louis  University  in 
the  early  1970s  using  readings  exclusively  from  analog  seismograms  first  focused  atten¬ 
tion  on  using  the  Lg  phase  for  seismic  source  size  estimation.  With  the  emergency  of 
widespread  digital  recordings  in  the  1980s,  the  focus  shifted  to  automatic  digital  process¬ 
ing,  using  in  particular  RMS  measurements  of  digitally  filtered  traces  (typically  0.6-3.0 
Hz)  in  a  fixed  time  window  (typically  2  minutes).  This  method  allowed  reliable  estimates 
to  be  made  even  at  veiy  low  SNR,  by  using  a  noise  compensation  procedure.  Thus,  NOR- 
SAR  and  Grafenbeig  array  measurements  of  Lg  waves  from  the  Shagan  River  area  were 
sufficiently  precise  to  allow  a  systematic  P-Lg  magnitude  bias  to  be  identified  between  the 
NE  and  SW  parts  of  that  site.  Available  yield  data  from  Soviet  sources  indicate  that  Lg 
magnitudes  show  the  best  consistency  with  yield  for  explosions  from  this  area.  Most  of  the 
evidence  for  the  **  stability”  of  Lg  magnitudes  is  indirect.  By  pairwise  comparison  of  RMS 
Lg  for  a  number  of  different  source-receiver  paths,  it  has  been  demonstrated  that  single¬ 
station  RMS  Lg  can  be  used  to  estimate  relative  magnitude  with  a  remaricably  small  scat¬ 
ter  (0.02-0-03  in  mi,  units  orthogonally).  This  was  first  shown  for  Semipalatinsk  explo¬ 
sions,  but  has  recently  been  confirmed  also  for  explosions  at  Novaya  Zemlya.  The  latter 
observation  is  particularly  interesting,  since  many  of  the  paths  exhibited  significant  ”Lg 
blockage”  effects.  Lg  waveforms  cannot  at  present  be  modelled  with  the  same  quality  of 
fit  between  synthetics  and  data  that  has  been  attained  with  other  phases,  and  many  aspects 
of  Lg  generation  and  propagation  characteristics  are  still  not  well  understood.  Among  top¬ 
ics  that  need  further  study  are  Lg  blockage  and  scattering  effects  caused  by  tectonic  heter¬ 
ogeneities  and  the  effects  of  topography,  near  source  geology  and  depth  of  burial  on  Lg 
excitation.  Much  challenging  wodc  therefore  remains  to  be  done  in  this  field. 
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Preface 


Under  contract  No.  F49620-C-89-()038,  NTNF/NORSAR  is  conducting  re.search  within  a 
wide  range  of  .subjects  relevant  to  seismic  monitoring.  The  emphasis  of  the  research  pro¬ 
gram  is  on  developing  and  assessing  methods  for  processing  of  data  recorded  by  networks 
of  small-aperture  arrays  and  3-component  stations,  forevents  both  at  regional  and  teleseis- 
mic  distances.  In  addition,  more  general  seismological  research  topics  are  addressed. 

Each  quarterly  technical  report  under  this  contract  presents  one  or  several  separate  investi¬ 
gations  addressing  specific  problems  within  the  scope  of  the  statement  of  work.  Summa¬ 
ries  of  the  re.search  efforts  within  the  program  as  a  whole  are  given  in  annual  technical 
reports. 

This  Scientific  Report  No.  1 1  is  the  annual  technical  report  for  the  period  1  October  1990  - 
30  September  1991.  It  contains  four  separate  contributions,  and  also  abstracts  for  the 
investigations  submitted  as  quarterly  technical  reports  during  FY91. 
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1  Summary 

This  Annual  Technical  Report  describes  the  work  accomplished  under  Contract  No. 
F49620-C-89-0038  during  the  period  1  October  1990  -  30  September  1991.  The  report 
contains  four  separate  contributions  (paragraphs  2. 1  through  2.4),  and  in  addition  ab.stracLs 
of  the  investigations  submitted  as  quarterly  technical  reports  during  FY91  (paragraph  2.5). 

The  inversion  of  NORSAR  travel  time  data  constituted  one  of  the  first  examples  of  seis¬ 
mic  tomography.  It  is  now  becoming  apparent  that  these  data  are  influenced  by  wave  dif¬ 
fraction  effects,  and  hence  that  ordinary  (ray  based)  tomography  may  be  inadequate.  In 
paragraph  2. 1  we  apply  a  newly  developed  tomographic  method  that  takes  diffraction 
effects  into  account.  Data  from  the  NORSAR  array  in  its  original  form  (22  subarrays; 
aperture  1(X)  km)  have  been  used  and  results  are  compared  with  those  of  the  ray-based 
inversion  experiments  undertaken  around  15  years  ago.  To  facilitate  comparison  with  ear¬ 
lier  results,  we  inverted  for  velocity  structure  in  the  same  box-like  structure  as  was  used  in 
the  previous  investigations.  As  we  Jo  not  know  the  true  structure  beneath  NORSAR,  we 
cannot  determine  the  model  fit  of  the  results,  but  the  data  fit  can  be  computed  by  forward 
modelling.  This  procedure  shows  that  the  travel  time  and  phase  delay  data  are  explained 
quite  well  by  both  methods.  On  the  other  hand,  ordinary  tomography  is  not  doing  well  in 
explaining  the  amplitudes,  whereas  diffraction  tomography  explains  about  20%  even 
when  amplitude  data  are  not  u.sed  in  the  inversion,  and  about  40%  when  the  inversion  does 
include  these  data.  It  is  concluded  that  the  modified  method  of  diffraction  tomography 
described  in  this  paragraph  promises  to  become  an  important  tool  for  determining  the 
structure  beneath  seismic  arrays. 

In  paragraph  2.2,  the  problem  of  discriminating  between  earthquakes  and  underground 
explosions  is  formulated  as  an  exercise  in  pattern  recognition  approach  analysis.  The  anal¬ 
ysis  has  been  applied  to  a  learning  set  of  44  nuclear  explosions  (8  test  sites)  and  35  earth¬ 
quakes  in  Eurasia  recorded  at  the  NORESS  array.  The  signal  features  considered  were  the 
normalized  power  in  8  spectral  bands  in  the  0.2-5.0  Hz  range  of  the  P-wave  (6  sec)  and  the 
P-coda  (30  sec).  Physically,  it  means  that  we  exploit  potential  differences  in  the  shape  of 
earthquake  and  explosion  spectra,  respectively.  Other  features  included  are  peak  P  and  P- 
coda  amplitude  frequencies  and  relative  P/P-coda  power.  These  19  features  were  extracted 
either  from  conventional  array  beam  traces  or  the  optimum  group  filtered  traces  (OGF- 
removal  of  coherent  low-frequency  noise).  Using  the  feature  selection  algorithm,  based  on 
estimates  of  the  expected  probability  of  misclassification  (EPMC),  only  2  to  4  features 
were  needed  for  optimum  discrimination  performance.  The  dominant  features  were  coda 
excitation  and  P-  and  P-coda  power  at  lower  signal  frequencies.  Furthermore,  feature 
parameters  extracted  from  the  OGF  traces  had  a  slightly  better  performance  in  comparison 
to  those  extracted  from  beam  traces.  Finally,  there  were  no  misclassifications  for  OGF- 
derived  features  when  the  explosion  population  was  limited  to  E.  Kazakh,  while  including 
events  from  the  other  test  sites  lead  to  a  decrease  in  discrimination  power. 

We  have  previously  applied  the  ray  tracing  technique  of  Bostock  and  Kennett  (1990)  to 
the  characterization  of  Lg  propagation  to  the  NORESS  array.  In  paragraph  2.3,  we  have 
considered  the  application  of  the  same  approach  of  representing  the  major  features  of  the 
Lg  phase  via  constructively  interfering  S  waves  multiply  reflected  within  the  earth’s  crust 
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to  understanding  azimuth  anomalies  observed  at  the  ACCESS  array.  The  motivation  for 
the  theoretical  modelling  has  been  the  rapid  changes  in  both  size  and  sign  of  ARCESS  Lg 
azimuth  anomalies  with  geographical  position,  clearly  observed  in  IMS  analysis  of  events 
to  the  south  and  east  of  Finland  (Gulf  of  Estonia  -  Leningrad  -  USSR  north  of  Lenihgrad). 
Theoretical  azimuth  anomalies  calculated  from  this  ray  tracing  scheme,  and  based  on  a 
simplified  model  of  the  Moho  structure  under  the  Fennoscandian  region,  are  found  to  be 
generally  consistent  with  the  trends  in  the  observed  anomalies  at  ARCESS.  However,  the 
ray  picture  for  Lg  is  probably  more  valuable  for  providing  a  tool  to  investigate  the  sensi¬ 
tivity  of  guided  wave  propagation  to  the  relative  position  of  the  source  and  major  subsur¬ 
face  features  .such  as  the  significant  Moho  topgraphy  under  Finiand. 

Paragraph  2.4  is  a  contribution  on  the  use  of  Lg  as  a  yield  estimation  tool.  In  recent  years, 
the  Lg  phase  has  emerged  as  maybe  the  most  promising  tool  to  obtain  precise  yield  esti¬ 
mates  by  seismic  means.  The  pioneering  work  by  Otto  Nuttli  of  Saint  Louis  University  in 
the  early  1970s  using  readings  exclusively  from  analog  seismograms  first  focused  atten¬ 
tion  on  using  the  Lg  phase  for  seismic  source  size  estimation.  With  the  emergency  of 
widespread  digital  recordings  in  the  1980s,  the  focus  shifted  to  automatic  digital  process¬ 
ing,  using  in  particular  RMS  measurements  of  digitally  filtered  traces  (typically  0.6-3.0 
Hz)  in  a  fixed  time  window  (typically  2  minutes).  This  method  allowed  reliable  estimates 
to  be  made  even  at  very  low  SNR,  by  using  a  noise  compensation  procedure.  Thus,  NOR- 
S  AR  and  Grafenberg  array  measurements  of  Lg  waves  from  the  Shagan  River  area  were 
sufficiently  preci.se  to  allow  a  sy.stematic  P-Lg  magnitude  bias  to  be  identified  between  the 
NE  and  SW  parts  of  that  site.  Available  yield  data  from  Soviet  sources  indicate  that  Lg 
magnitudes  show  the  best  consistency  with  yield  for  explosions  from  this  area.  Mo.st  of  the 
evidence  for  the  “stability”  of  Lg  magnitudes  is  indirect.  By  pairwise  comparison  of  RMS 
Lg  for  a  number  of  different  source-receiver  paths,  i*.  has  been  demonstrated  that  single¬ 
station  RMS  Lg  can  be  used  to  estimate  relative  magnitude  with  a  remarkably  small  scat¬ 
ter  (0.02-0-03  in  m^  units  orthogonally),  litis  was  first  shown  for  Semipalatin.sk  explo¬ 
sions,  but  has  recently  been  confirmed  also  for  explosions  at  Novaya  Zemlya.  The  latter 
ob.servation  is  particularly  intere.sting,  since  many  of  the  paths  exhibited  significant  “Lg 
blockage”  effects.  Lg  waveforms  cannot  at  pre.sent  be  modelled  with  the  same  quality  of 
fit  between  synthetics  and  data  that  has  been  attained  with  other  phases,  and  many  aspects 
of  Lg  generation  and  propagation  characteristics  are  still  not  well  understood.  Among  top¬ 
ics  that  need  further  study  are  Lg  blockage  and  scattering  effects  caused  by  tectonic  heter¬ 
ogeneities  and  the  effects  of  topography,  near  source  geology  and  depth  of  burial  on  Lg 
excitation.  Much  challenging  work  therefore  remains  to  be  done  in  this  field. 
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2  Summary  of  Technical  Findings  and  Accomplishments 

2.1  Seismic  diffraction  tomography  of  array  data 

Introduction 

Seismic  tomography  has  become  a  widely  applied  method  to  determine  the  crustal  and 
lithospheric  velocity  structure  beneath  an  array.  Actually,  different  inversion  algorithms 
are  being  applied,  but  they  all  invert  travel  time  data  by  linearizing  the  problem,  appealing 
to  Fermat’s  principle.  Apart  from  its  inherent  scientific  interest,  the  result  can  be  used  to 
interpolate  travel  time  residuals;  as  such,  tomographic  results  are  also  of  operational  inter¬ 
est. 

Despite  the  widespread  application  of  the  method,  a  number  of  limitations  were  recog¬ 
nized  at  an  early  stage.  A  fundamental  limitation  is  due  to  the  neglect  of  wave  diffraction 
off  the  geometrical  ray  path.  This  neglect  limits  the  resolution  of  tomographic  results  in 
general.  The  limitation  is  especially  serious  in  circumstances  where  diffraction  effects 
dominate  the  observations,  as  may  be  expected  in  the  presence  of  low  velocity  regions 
among  others.  Another  practical  problem  is  that  seismic  tomography  operates  under  the 
a.ssumption  of  travel  time  data  representing  first  arrivals,  but  these  are  usually  not  observ¬ 
able.  Instead  one  measures  detection  times,  or  alternatively  at  an  array,  one  may  identify 
relative  time  delays  with  the  waveform  correlation  lags  between  sensors.  In  either  case, 
diffraction  may  significantly  affect  the  data.  Diffraction  effects  have  been  taken  into 
account,  in  a  first-order  approximation,  in  methods  called  diffraction  tomography  (e.g., 
Devaney,  1984).  Unfortunately,  the  first-order  approximation  does  not  adequately  recover 
the  phase  and  hence  the  travel  time  residual  (although  the  so-called  Rytov  approximation 
has  been  designed  to  improve  this).  Another  general  problem  is  that  non-physical  diffrac¬ 
tion  from  the  boundary  of  the  region  under  study  may  interfere  with  the  solution  (unless 
the  velocity  perturbation  is  zero  outside  this  region).  Finally,  current  methods  of  diffrac¬ 
tion  tomography  assume  a  homogeneous  background  medium  and  regular  sampling  of  the 
wave  field,  such  that  Fourier  transform  methods  can  be  used  for  the  solution  of  the  scatter¬ 
ing  problem  at  hand.  In  its  prese.>it  form  these  methods  are  of  limited  use,  at  least  in  earth¬ 
quake  seismology. 

We  have  experimented  with  a  new  formulation  that  combines  diffraction  with  ordinaiy 
seismic  tomography  (Doombos,  1991),  overcoming  the  problems  mentioned  above.  We 
anticipate  that  the  method  can  be  applied  to  a  range  of  problems,  including  for  example 
cross-well  tomography.  However,  the  experiment  we  describe  here  involves  data  from  the 
NORSAR  array  in  its  original  (laige  aperture)  configuration.  It  is  of  interest  to  note  that 
the  travel  time  data  from  this  array  have  been  used  in  one  of  the  first  applications  of  seis¬ 
mic  tomography  (Aki  et  al,  1977),  although  at  the  time  the  inverse  method  was  not  called 
tomography.  Thus  our  experiment  naturally  leads  to  a  comparative  study  of  diffraction  and 
ordinary  tomography  applied  to  teleseismic  array  data. 
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Outline  of  theoretical  background 

The  basis  of  ordinary  seismic  tomography  is  the  simple  relation,  in  the  ray  approximation: 

8t=  I  5^  da  (1) 

ray 

Here  5t  is  the  change  of  travel  time  t  due  to  a  slowness  perturbation  8s  along  the  ray,  and 
using  Fermat’s  principle,  the  integration  is  along  the  ray  path  in  the  unperturbed  reference 
structure.  The  ray  approximation  thus  predicts  a  time  shift  of  the  reference  pulse: 

m(/)  =m^(/-8t)  (2) 

Diffraction  theory,  on  the  other  hand,  predicts 

M(r)  =M^(r) +8m(/)  (3) 

and  the  linear  (Bom)  approximation  for  8m  is  of  the  form 

8m  (/)  =  I  F  (c,  p,  G,  M„)  8s  dV  (4) 

V 

For  example,  for  an  acoustic  medium  with  velocity  c.  density  p  and  Green’s  tensor  Q  : 

F  =  -2pcG*M^  (5) 

Upon  inspecting  equations  (2),  (3)  and  (4).  one  may  anticipate  the  following  problems: 

(1 )  If  the  perturbed  structure  is  relatively  smooth,  permitting  the  ray  approximation  (2), 
then  from  (3): 


8m  (r)  =-m^(/)8t 

and  this  is  obviously  a  poor  approximation  if  8r  is  no  longer  small. 

(2)  Inverting  equation  (4)  implies  the  assumption  8s  =  0  outside  V,  causing  numerical 
diffraction  if  in  fact  8s  ^  0  outside  V. 

In  view  of  the  above  it  certainly  would  be  an  advantage  to  relate  the  perturbation  term  to 
the  time-shifted  pulse: 


yit)  =  M^(t-5T) +8M^(r-8i)  (6) 

provided  there  exists  a  sufficiently  simple  relation  between  8m^  and  8s.  Such  a  relation 
was  derived  by  Doombos  (1991 ): 


•  V8s 


\dV 


(7) 
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where  p°,p^  are  the  slowness  vectors  of  the  incident  and  scattered  wave,  respectively. 
Equation  (7)  is  linear  in  the  slowness  gradient,  and  any  of  the  velocity  or  slowness  param¬ 
eterization  schemes  in  common  use  will  convert  it  to  a  form  that  is  linear  in  the  slowness 
itself,  so  that  standard  inverse  methods  can  be  applied.  The  common  case  of  a  block 
parameterization  is  special  in  the  sense  that  the  diffraction  term  contributes  only  diffrac¬ 
tions  from  the  block  boundaries: 


(8) 


where  is  the  boundary  between  block;  and  k  with  normal  n,/^,  the  sum  over  k  is 
restricted  to  blocks  adjacent  to  j,  and  the  additional  restriction  K>j\&  used  to  avoid  dupli¬ 
cating  interfaces. 

Writing  the  perturbation  in  the  form  (6)  is  consistent  with  the  Rytov  approximation.  The 
approximation  was  derived  in  the  frequency  domain: 


In 


t/((D) 


=  /(i)5t  + 


8Uj{(P) 


(9) 


where  8t  is  given  by  equation  (1),  and  by  the  Fourier  transform  of  equation  (7)  or 
(8).  Equation  (9)  together  with  (8)  is  linear  in  the  slowness  perturbation  5s;  in  practice  a 
solution  must  be  obtained  by  iteration.  Several  iteration  schemes  are  possible,  but  in  any 
case  it  is  important  to  up.  date  the  phase  delay  of  the  diffracted  waveheld  in  (7).  In  the 
present  work  we  use 


In 


s  /(i)5t  + 


5f/^(0)) 


(10) 


where  Uj  ((u)  is  the  diffracted  field  after  the  previous  iteration,  and  the  phase  delay  in 
and  is  updated  using  the  approximation  (1). 


A  model  experiment 

The  model  consists  of  a  block-like  structure  of  velocity  perturbations.  The  block  size  is  20 
km,  and  the  velocity  varies  continuously  by  ±  3%  between  blocks,  over  a  gradient  zone 
of  about  6  km  (Fig.  2.1.1).  The  size  of  the  whole  structure  is  2(K)  x  200  km  horizontally, 
and  40  km  vertically.  The  top  of  the  structure  is  at  80  km  depth.  An  array  with  the  config¬ 
uration  of  NORS  AR  is  assumed  to  record  tele.seismic  P  waves  coming  from  different 
directions,  as  shown  by  the  slowness  diagram  of  Fig.  2.1.2.  The  total  data  set  is  given  by 
the  P  waves  from  40  events  recorded  at  22  stations.  The  incident  pulse  form  is  chosen  to 
be  similar  to  a  typical  teleseismic  record.  Synthetic  waveforms  were  generated  by  a  com¬ 
bination  of  3-D  dynamic  ray  tracing  and  Kirchhoff-type  integration:  Rays  were  traced 
through  the  3-D  structure  generating  the  wavefield  on  an  integration  surface  just  above  the 
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top  of  the  Structure,  The  field  at  the  surface  was  then  determined  by  a  Kirchhoff-lype  of 
integral  over  this  surface.  Relative  travel  time  residuals  were  determined  by  correlating 
the  waveforms  at  the  22  stations  with  an  array  beam.  These  residuals  form  the  data  set  for 
ordinary  tomography  based  on  equation  (1),  The  waveform  spectrum  f/  (to)  was  also 
determined  for  each  station,  together  with  the  array  beam  spectrum  (w) ,  The  imagi* 
n{U7  part  of  In  ( U/U^)  in  equation  (9)  gives  the  phase  delay  spectrum,  which  is  the  coun¬ 
terpart  of  the  travel  time  residual  5t,  Inversion  of  the  phase  delay  spectra  (that  is, 
operating  on  the  imaginaiy  part  of  equation  (9))  constitutes  the  central  pan  of  diffraction 
tomography.  In  this  experiment  we  inverted  for  6  frequencies  in  the  band  0,6-2, 1  Hz, 
Simultaneous  inversion  of  the  amplitude  spectra  (the  real  part  of  equation  (9))  Is  expected 
to  increase  stability  of  the  results. 

Obviously,  the  model  considered  here  is  not  meant  to  be  a  realistic  model  of  Earth  struc¬ 
ture,  Rather  it  serves  as  a  test  to  evaluate  the  inversion  algorithm,  and  it  allows  a  compari¬ 
son  of  different  tomographic  methods  under  controlled  circumstances.  The  model  is 
special  in  its  regular  sequence  of  low  and  high  velocity  regions.  Diffraction  around  the 
low  velocity  regions  tends  to  mask  the  signal  from  such  regions  in  the  travel  time  data. 
This  process  is  particularly  effective  if  the  size  of  a  region  is  of  the  order  of  a  Fresnisl 
zone;  the  scaling  of  our  experiment  was  in  fact  designed  to  match  this  condition.  Another 
special  feature  of  our  experiment  is  that  in  the  inversion  we  choose  the  block  parameter¬ 
ization  to  match  the  block  size  structure  of  the  model.  This  will  minimize  the  smoothing  of 
tomographic  results,  but  it  will  also  emphasize  errors  in  these  results. 

Under  these  special  conditions,  ordinary  tomography  does  indeed  produce  relatively  large 
errors.  We  define  a  measure  of  the  normalized  squared  error: 

e  =  J  (5c-6c)  W/J  (8c)2</V'  (ID 

V  V 

where  8c  is  the  model  velocity  perturbation  and  Sc  is  its  estimate.  The  values  for  £ 
depends  on  how  strongly  we  damp  the  solution,  but  for  ordinary  tomography  we  always 
find  e  >  1(K)% .  For  diffraction  tomography  we  get  e  =  35% .  The  error  is  largest  near  the 
boundary  of  the  model.  Deleting  from  integration  (1 1)  the  blocks  along  the  boundary 
reduces  £  to  25%.  This  represents  a  significant  improvement  over  ordinary  tomography. 
For  other,  smoother  models,  the  improvement  would  of  course  be  less  dramatic. 

Application  to  the  NORSAR  array 

The  geometry  of  the  experiment  is  illustrated  in  Fig.  2. 1 .3.  To  facilitate  comparison  to  ear¬ 
lier  results,  we  invert  for  velocity  suncture  in  the  same  box-like  region  as  Aki  et  al  (1977) 
and  we  use  the  same  block  parameterization.  However,  to  more  fully  include  the  diffracted 
field  we  extended  the  box  size  slightly  horizontally.  The  size  was  220  x  220  km  horizon¬ 
tally  and  126  km  vertically,  layer  thicknesses  were  17, 19, 30, 30  and  30  km.  and  the  block 
size  within  each  layer  was  20  x  20  km  horizontally:  We  also  use  the  same  basic  data  set, 
retainin j  the  data  from  86  events.  The  P  slowness  solutions  from  these  events  are  shown 
in  Fig.  1.1.4,  Clearly,  sampling  of  the  structure  is  very  non-uniform;  in  particular  the 
southern  part  of  the  structure  is  poorly  sampled  by  rays.  It  is  well-known  that  non-uniform 
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sampling  tends  to  deteriorate  tomographic  results,  but  this  problem  would  be  less  severe 
for  diffraction  tomography  due  to  the  additional  constraint  by  the  diffraction  term  (equa¬ 
tion  7).  An  example  of  the  result  for  one  of  the  layers  (the  depth  range  66-96  km)  is  shown 
in  Fig.  2.1.5  for  ordinary  tomography,  and  in  Fig.  2.1.6  for  diffraction  tomography.  Ignor¬ 
ing  the  prominent  features  near  the  boundary  of  the  grid  as  possibly  due  to  numerical 
boundary  diffraction,  the  pattern  near  the  center  of  the  grid  shows  about  ±  6%  velocity 
perturbations  for  both  methods.  Details  of  the  two  results  are  different  however.  We  do  not 
know  the  true  structure  so  we  cannot  determine  the  model  fit  of  the  results,  but  we  can 
compute  the  data  fit.  Synthetic  data  were  constructed  after  forward  modelling  the  dif¬ 
fracted  field  using  a  first-order  (Rytov)  approximation  and  adjusting  the  pha.se  delay.  A 
summary  of  results  is  shown  in  Fig.  2. 1 .7a  for  different  types  of  data  (travel  time  residuals, 
phase  delays  and  amplitudes).  Clearly  the  travel  time  and  phase  delay  data  arc  explained 
quite  well  by  both  methods.  On  the  other  hand,  ordinary  tomography  is  not  doing  well  in 
explaining  the  amplitudes;  diffraction  tomography  explains  about  20%  even  when  ampli¬ 
tude  data  are  not  u.sed  in  the  inversion,  and  about  40%  when  the  inversion  does  include 
the.se  data.  In  a  further  test  we  inverted  the  synthetic  data  for  one  of  the  models  (from  dif¬ 
fraction  tomography:  DT,  p&a  of  Fig.  2.1.7a).  This  produces  a  measure  of  model  fit  in 
addition  to  the  data  fit,  for  the  different  methods  of  tomography.  The  results  are  summa¬ 
rized  in  Fig.  2. 1 .7b.  The  general  features  of  the  data  fit  are  similar  to  those  for  real  data. 
Ordinary  tomography  produces  a  very  poor  fit  to  the  true  model,  but  diffraction  tomogra¬ 
phy  reduces  the  model  fit  error  e  (equation  1 1 )  to  about  30%.  Of  course  this  may  be  an 
overly  optimistic  result  since  the  forward  and  inverse  method  were  based  on  the  same 
approximation,  but  the  results  do  agree  with  those  of  the  model  experiment  described  ear¬ 
lier.  The  modified  method  of  diffraction  tomography  described  here  thus  promises  to 
become  an  imponant  tool  for  determining  the  structure  beneath  seismic  arrays. 


E.  0degaard 
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Hg.  2.1.1.  Model  velocity  structure,  with  location  of  receiver  array  indicated.  The  contin¬ 
uous  velocity  viiriation  between  blocks  is  not  shown  in  this  figure. 
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Fig.  2.13.  Sketch  of  inversion  experiment  at  NORSAR. 
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Sloymess  solutions  for  the  86  events 
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Fig.  2.1.4.  Slowness  of  P  waves  used  for  inversion  experiment  at  NORSAR. 


11 


Basic  Scismological  Research 


November  1991 


Fig.  2.1.5.  Velocity  perturbations  for  layer  4  (66-96  km).  Results  from  onlinary  tomogra¬ 
phy.  Open  fields  are  not  sampled  by  rays. 
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Relative  velocity  perturbation 
Layer  4  (SS—QSkn) 


Fig.  2.1.6.  Velocity  penurbations  for  layer  4  (66-96  km).  Results  from  diffraction  tomog¬ 
raphy. 
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Fig.  2.1.7a  (left):  Data  fit  after  tomography  at  NORSAR.  ST:  Ordinary  tomography.  DT,  p 
and  DT,  p&a:  Diffraction  tomography  without  and  with  amplitude  data,  b  (right):  Data 
and  model  fit  after  tomography  with  synthetic  data  at  NORSAR.  Synthetics  were  gener¬ 
ated  for  model  DT,  p&a. 
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2.2  Enhanced  seismic  source  discrimination  using  NORESS 
recordings  from  Eurasian  events 


Introduction 

The  problem  of  distinguishing  underground  nuclear  explosions  from  natural  earthquakes 
using  seismic  data  has  been  studied  for  a  long  time.  In  fact,  it  dates  back  to  19  sep  1957  as 
on  that  day  a  nuclear  explosion,  code-named  Rainier,  was  detonated  under  the  Nevada 
desert.  In  order  to  avoid  the  once  troublesome  issue  of  in-country  operation  of  non¬ 
national  seismic  stations,  the  source  identification  research  in  the  1960  and  1970-ties 
focused  on  observations  in  the  tcleseismic  window.  The  most  successful  criteria  for  seis¬ 
mic  source  identification  were  spectral  ratio  variants,  mainly  between  non-overlapping 
frequency  bands  for  the  P-signal  itself,  or  the  relative  signal  excitation  at  1  sec  (P-wave) 
and  20  sec  (surface  waves).  The  latter  is  often  denoted  the  mb'.Mg  (body  wave  versus  sur¬ 
face  wave  magnitudes)  discriminant.  Another  commonly  used  discriminant  was  the  so- 
called  complexity  tied  to  the  ratio  of  P  coda  RMS  in  two  consecutive  windows  of  lengths 
around  5  sec  and  30  .sec  respectively  (Dahlman  and  Israelson,  1977).  A  variant  of  the  com¬ 
plexity  often  denoted  the  coda  discriminant,  was  introduced  by  Tjostheim  (1975, 1978) 
and  Sandvin  and  Tj0stheim  (1978)  using  autoregressive  spectral  coefficients  as  signal 
attributes  in  combination  with  more  advanced  discrimination  statistics.  Recently,  new  con¬ 
cepts  have  been  introduced,  namely,  the  so-called  artificial  or  trained  neural  networks 
technique  (e.g.  see  Dowla  et  al,  1990  and  Dysart  and  Pulli,  1990). 

In  this  study  we  address  the  problem  of  teleseismic  source  discrimination  and  explore  the 
potential  of  the  spectral  ratio  and  complexity  discriminants.  We  u.se  data  from  the  NOR¬ 
ESS  array,  which  has  an  excellent  detectability  for  events  in  parts  of  Eurasia. 

Event  selection  and  NORESS  record  preprocessing 

The  pre.sumed  earthquakes  (PDE  and  ISC  listings)  and  presumed  underground  nuclear 
explosions  (NORSAR  listings)  used  in  this  .study  are  listed  in  Tables  2.2.1  and  2.2.2  and 
depicted  in  Fig.  2.2.1  as  well.  Most  of  the  explosions  (32  out  of  44)  stem  from  the  E. 
Kazakh  (Semipalatinsk)  test  site  and  hence  the  earthquakes  in  the  m^  range  4.0  to  6.0  were 
selected  from  the  same  general  area.  The  other  test  sites  are  in  aseismic  areas  for  which 
sufficient  earthquake  recordings  are  not  available.  However,  the  outlaying  explosions  were 
included  for  a  check  on  discriminant  robustness.  A  simple  and  efficient  scheme  for  SNR 
enhancement  is  delay-and-sum  processing  or  beamforming  which  is  most  efficient  in  the 
2-8  Hz  band.  At  lower  frequencies  (below  say  2  Hz)  the  wavelengths  of  microseisms  are 
of  the  same  order  as  the  array  aperture,  and  hence  strong  correlation  in  the  noise  across  the 
a* ray  is  often  observed.  In  .such  cases  Ingate  et  al  (1985)  demonstrated  the  usefulness  of 
niHximum,  likelihood  schemes  for  suppressing  correlated  noise.  Recently,  even  more 
iidvanc^J  methods  have  been  inU'oduced  by  Kushnir  et  al  (1990),  which  are  extensively 
used  here  for  suppressing  low  frequency  propagating  noise.  The  efficiency  of  the  this 
scheme,  commonly  denoted  the  Optimal  Group  Filtering  (OGF)  technique  is  demon¬ 
strated  in  Fig.  2.2.2.  Naturally,  removal  of  low  frequency  noise  is  important  as  both  theo¬ 
retical  and  observational  studies  demonstrate  that  part  of  the  discrimination  power  is 
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vested  in  the  low/high  fre(|uency  bands  of  the  P-signal  (Evemden  et  al,  1986).  Evemden 
(1977)  advocates  the  use  of  P-signal  frequencies  up  to  9  Hz  in  teleseismic  discrimination 
studies,  while  we  restricted  the  analysis  to  5  Hz.  Our  rationale  being  that  there  is  not  much 
signal  energy  above  5  Hz  and  besides  cultural  sources  like  local  quarry  blasting,  mining 
explosions  and  fast  running  machinery  (saw  mills)  could  easily  bias  the  observational 
data.  In  short,  our  discriminant  parameters  were  extracted  both  from  single  P-beam  traces 
and  from  OOF  traces  in  order  to  have  observations  forjudging  the  relative  importance  of 
the  latter. 

Signal  attributes  -  class  of  wavefield  parameters  for  discrimination 

As  mentioned,  the  most  powerful  discrimination  parameters  are  related  to  differential  sig¬ 
nal  excitation  in  different  frequency  bands  for  earthquakes  and  explosions  respectively.  P- 
wave  parameters  are  an  obvious  choice  here,  because  this  phase  is  most  easily  detected. 
Coda  waves  are  interesting  as  they  not  only  reflect  the  source  but  also  the  source  location 
within  the  crust/upper  mantle.  For  example,  coda  excitation  and  duration  are  far  less  effi¬ 
cient  for  surface  explosions  and  deep  earthquakes  relative  to  shallow  and  intermediate 
depth  earthquakes  (Dainty.  1990).  Rayleigh  waves  are  not  considered  simply  because  they 
become  embedded  in  background  noise  for  event  magnitudes  at  4.5  -  5.0  or  below.  The 
essence  of  the  above  discussion  is  that  our  class  of  discrimination  ptirameters  is  tied  to  the 
P-signal  and  its  coda  as  illustrated  in  >^igs.  2.2.3  and  2.2.4.  The  choice  of  window  lengths 
of  6  .sec  and  30  sec  reflects  that  most  of  the  desirable  information  from  teleseismic.  but 
also  regional  events  tire  accumulated  in  the  first  3  -  6  sec  of  the  P-wave  and  then  the  sub.se- 
(|ucnt  10  -  30  sec  coda  waves  where  scattering  contributions  are  still  significant. 

Proce.ssing  details  were  as  follows.  For  each  event  5  minutes  of  recordings  for  all  25  verti¬ 
cal  NORESS  sensors  were  extracted.  The  first  2  minutes  of  pure  noise  were  used  for  esti¬ 
mating  OGF-filter  coefficients  for  the  slowne.ss  and  azimuth  of  the  individual  events.  After 
the  filtering  was  performed,  amplitude  spectra  (FFT)  were  calculated  for  the  P-signal  (6 
.sec)  and  the  P  coda  (30  sec).  The  power  spectrum  for  the  non-overlapping  8  bands  speci¬ 
fied  in  Figs.  2.2.3  were  obtained  by  simple  averaging  of  spectral  squared  amplitudes. 
Since  the  events  used  in  analysis  have  widely  different  magnitudes,  all  power  spectfa  were 
normalized  by  their  maxinuims.  Physically,  this  means  that  potential  spectral  shape  differ¬ 
ences  arc  exploited  for  event  discrimination  purposes.  The  final  feature  parameter  values 
were  obtained  by  taking  the  logarithm  of  the  normalized  spectral  values.  As  shown  in 
Fig.2.2.3  we  have  8  feature  parameters  for  both  the  P  signal  and  the  Pcoda  for  the  same 
set  of  frequency  bands.  Additional  panuneters  introduced  (nos.  17  and  18)  are  peak  spec¬ 
tral  frequencies  and  finally  the  19th  parameter  being  the  ratio  of  P/P  coda  .speco^l  max¬ 
ima.  Note  that  the  last  parameter  is  close  to  the  classical  complexity  definition. 

Di scr i mi na t ion  approach 

The  seismic  discrimination  undertaking  is  a  two-stage  process;  firstly,  relevant  discrimina¬ 
tion  parameters  must  be  defined  and  extracted  from  the  records,  and  secondly,  decision 
rules  (discriminators)  must  be  introduced  to  ensure  proper  event  classification.  At  the  out¬ 
set  of  this  study,  we  discussed  rather  extensively  among  ourselves  which  discrimination 
approach  would  be  best  for  classification  of  Eurasian  events  as  recorded  by  NORESS.  The 
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importance  of  the  physical  aspect  of  the  problem  at  hand  was  duly  recognized,  that  is,  the 
extracted  discrimination  parameters  must  be  seismologically  relevant.  Furthermore,  the 
persistence  of  the  seismic  source  identification  problem  is  that  a  unique  set  of  discrimi¬ 
nants  for  weak  events  apparently  does  not  exist,  so  we  wanted  to  have  the  ability  to  decide 
statistically  which  ones  of  such  parameters  would  be  most  informative  for  a  given  area  or 
region. 

The  recognition  problem  in  our  case  may  be  formulated  as  follows.  The  data  set  consists 
of  three  sets.  The  first  two  are  the  sets  R,  of  explosions  (NE)  spectral  parameters  vectors 
and  R2  earthquake  (EQ)  spectral  parameters  vectors  r/^  : 


{r 


O') 


fO),  . 


rj-'^  =  (p,,...,p„);  /=1,2 


(I) 


Nj  -  number  of  observations  in  the  j-th  learning  set  (earthquakes  or  explosions),7=7,2 
Pf.-  values  of  classification  features  (spectral  parameters  P  and  P-coda),  k=:l,n 


n  -  number  of  classification  features. 


R,  and  Rj  form  the  learning  material  for  the  NE  and  EQ  classification  features.  The  third 
.set  is  the  vector  X  containing  discrimination  parameters  for  an  “unknown”  event,  i.e.,  the 
vector  being  classified  with  no  prior  knowledge  as  to  its  .source  identification.  On  the  basis 
of  the  learning  samples  R,  and  R2.  we  must  make  a  decision  on  whether  X  belongs  to  the 
first  or  the  second  class.  This  is  done  by  using  a  decision  (discriminator)  function 
g  (X,  Rp  R2)  •  The  equivalent  Bayesian  rules  are:  If  g  (X,  R,,  R2)  <  C,  the  hypothesis 
HI  is  correct  and  X  belongs  to  the  first  class;  if  g  (X,  R,,  R2)  >  C,  the  hypothesis  H2  is 
correct  and  X  belongs  to  the  .second  class.  Here  C  is  a  constant,  that  is,  a  threshold  value. 
Using  this  decision  rule,  there  are  errors  of  two  kinds;  the  vector  X  belongs  to  the  second 
class  (the  hypothesis  H2  is  correct),  while  it  is  assigned  to  the  first  class.  Errors  of  the  sec¬ 
ond  kind  represent  the  reverse  situation.  Errors  occur  because  the  observations  Rj  ,R2  and 
X  are  random,  measurement  errors  and  the  random  nature  of  the  discrimination  features 
themselves.  The  classification  errors  of  the  first  and  second  kind  are  described  by  the  prob¬ 
abilities  PI  and  P2.  It  is  also  assumed  that  a  priori  probabilities  of  the  vector  X 

belonging  to  the  j-th  class  jire  known,  and  the  total  error  probability  +  ^2^2 

can  then  be  evaluated.  In  our  case  the  sizes  of  the  learning  samples  (3S  EQ  and  44  NE)  are 
comparable  to  the  dimension  n  on  of  the  features  space  ( 19  spectral  parameters).  In  such 
cases,  it  is  difficult  to  resolve  the  basic  classification  problem  -  to  select  the  optimum 
decision  rule  and  evaluate  its  corresponding  error  probabilities.  Since  the  distribution  pat¬ 
terns  corresponding  to  the  first  and  second  class  are  unknown,  it  is  theoretically  impossi¬ 
ble  to  construct  a  uniformly  optimum  decision  rule,  which  in  all  cases  will  yield  the  least 
probability  of  misclassification  (Kushnir  et  al,  1986).  A  criterion  for  choosing  a  par¬ 
ticular  rule  is  the  requirement  of  minimizing  the  probability  of  misclassification  (PMC). 

From  this  point  of  view,  when  the  Rj  and  R2  di.stributions  are  normal,  the  linear  di.scrim- 
inant  function  (LDF)  is  optimal  in  some  statistical  sense  (Anderson,  1958): 

/.(X.Rj.Rj)  =  (X-0.5(/-,+r2))*S''(/-,-r2)  <C  (2) 
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where 


X  -  vector  being  classified, 

S  *  covariance  matrix  of  learning  sets, 
fj  -  average  vector  of  j-th  learning  set, 

C  -  threshold;  above  is  EQ,  below  is  NE 

•  number  of  observations  in  the  j-th  learning  set 
*  -  matrix  transpose. 

Note,  the  LDF  expression  in  eq.  (2)  is  often  denoted  the  Fischer  discriminant  (Anderson, 
1958). 

Error  probability  estimation 

The  principal  method  of  deriving  approximate  expressions  for  error  probabilities  is  that  of 
asymptotic  representations  of  the  distributions  of  discriminator  statistics.  Of  major  interest 
in  seismic  discrimination  analysis  is  the  case  when  N  (number  of  events  in  learning  sets) 
and  n  (number  of  features)  are  of  the  same  order  and  have  double  digit  values.  For  LDF 
the  asymptotic  equation  suggested  by  Kolmogorov,  where  the  error  probabilities  are 
investigated  when  both  n  and  N  approach  infinity  and  n/N  is  constant,  remains  approxi¬ 
mately  valid.  Deev  (1970)  has  shown  that  in  Kolmogorov’s  asymptotic,  the  distribution  of 
the  LDF  linear  discriminator  (L  in  eq.  (2))  is  asymptotically  normal  with  mean  M  and  vari¬ 
ance  V: 


[2{N-l)/{2N-\~n)\[{-\)iV^/2]  j=  1,2 

V=  [{N-l){2N-\){2N+l)/{2N-n-l){2n-n)n][D^  +  2n/N]  (3) 

with  N  =  Ni=  N2,  and  =  (r,  -  r2)*  S’^ri  -  r2). 
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From  (3),  the  following  approximate  expression  for  the  EPMC  of  the  LDF  (Linear  Dis¬ 
criminant 


Function)  can  be  written: 


where 


PI  =P{j?<q//2}  =G((C-M,)/V) 
P2  =  P{i'>q//1}  =  G((A/2-C)/V) 

y 

G(y)  =  ]/j2n  J 


(4) 


Raudis  and  Pikyalis  (1975)  have  shown  that  the  relative  error  in  evaluating  the  probability 
of  misclassification  yielded  by  asymptotic  expression  (4),  as  compared  with  the  error  com¬ 
puted  from  an  exact  expression,  is  not  greater  then  few  percent  in  the  common  observa¬ 
tions  ranges  of  N  and  n.  The  event  classification  strategy  as  adopted  in  this  study  is 
visualized  in  Fig.  2.2.5. 

Broad  practical  applications  of  LDF  in  various  classification  problems  have  confirmed  its 
efficiency  for  small  sizes  N1  and  N2  of  learning  samples  and  for  distributions  that  differ 
from  normal  (Azen  et  al,  1975;  Weber  et  al.  1986;  Tsvang  et  al,  1986).  In  essence,  the 
advantage  of  LDF  is  its  simplicity,  and  the  exi.stence  of  a  method  for  calculating  EMPC. 

Selection  of  most  informative  features 

As  mentioned  above,  the  linear  di.scriminant  function  has  a  good  performance  when  the 
sizes  of  learning  sets  are  small.  LDF  can  be  applied  not  only  to  the  total  number  of  dis¬ 
criminant  features  (rj^^  in  eq.  (1)),  but  also  to  vectors  made  up  of  a  subset  of  the  discrimi¬ 
nants.  Various  values  of  the  expected  probability  of  misclassification  (eq.  (4)). 

The  natural  strategy  would  be  to  take  only  the  set  of  features  which  represents  the  mini¬ 
mum  of  the  EPMC.  In  order  to  make  use  of  this  idea,  one  must  be  able  to  arrange  the 
available  features  in  a  sequence  such  that  inaeasing  the  set  of  features  by  one  feature  at  a 
time  would  produce  the  maximum  rate  of  growth  of  the  function  D^(n).  One  possible  way 
of  ranking  the  features  is  the  following:  at  each  iteration  step  K  the  optimum  subset  of  K-J 
features  previously  selected  is  increased  by  adding  the  one  feature  from  among  the 
remaining  set,  which  yields  the  laigest  difference  d\K)~d\K-1).  The  features  are  thus 
ranked  in  significance  and  in  relation  to  the  Mahalonobis  distance  D  as  a  function  of  fea¬ 
tures  n.  Then  the  EPMC  are  computed  and  plotted  versus  ranked  features.  The  optimal 
subset  of  features  corresponds  to  the  minimum  of  EPMC,  as  illustrated  in  Fig.  2.2.6. 

To  provide  more  confident  decision  we  used  also  the  so-called  Jack-knife  method,  which 
allows  the  creation  of  an  independent  data  set  for  testing  the  discriminator  without 
decreasing  the  learning  data  set.  Here  learning  and  testing  are  carried  out  repeatedly:  dur¬ 
ing  each  learning  sequence  one  element  is  omitted  from  the  entire  sample  and  in  turn  used 
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for  testing.  The  total  number  of  errors  is  then  summarized.  In  essence,  all  elements  of  the 
sample  are  used  both  for  learning  and  testing,  but  the  element  being  used  for  recognition 
m'e  statistically  independent  of  the  learning  set.  Finally,  it  should  be  noted  that  the  Jack' 
Knife  method  may  also  be  used  for  estimating  the  conditional  probability  of  misclassifica- 
tion  not  the  EPMC. 

Results 

Perhaps  the  most  important  aspect  of  the  seismic  source  identification  problem  is  that  of 
finding  proper  discriminant  parameters.  In  our  case  we  selected  19  discriminant  candi¬ 
dates  (Fig.  2.2.3)  and  have  used  the  Expected  Probability  of  Misclassification  Probability 
(EPMC)  measure  for  ranking  the  relative  importance  of  these  parameters  as  illustrated  in 
Fig$.2.2.7a  and  2.2.7b.  In  these  figures  we  differentiate  between  two  populations,  namely, 
all  events  (44  NE  and  3.5  EQ)  and  only  Semipalatinsk  explosions  and  all  earthquakes  (32 
NE  35  EQ)  for  OGF  (Optimum  Group  Filtering)  and  BEAM  traces,  respectively.  Param¬ 
eter  no.  19  (Fig.  2.2.3)  are  clearly  the  most  dominant  for  all  cases,  and  thus  demonstrate 
the  ieje\"?nce  of  the  classical  complexity  discriminant.  There  appear  to  be  some  differ¬ 
ences  between  the  OGF  and  BEAM  trace  derived  parameters,  as  in  the  former  case  the 
optimum  performance  (EPMC  minimum)  is  tied  exclusively  to  the  coda  parameters.  For 
BEAM  traces,  the  classical  spectral  parameters  are  of  some  importance  as  weight  is  given 
to  relative  signal  power  in  the  low  frequency  bands  (parameters  1  and  3  in  Fig.  2.2.7a). 

The  parameters  no  10  and  19  derived  from  OGF  traces  have  clearly  the  best  performance 
without  any  misclassification  when  the  NE  populations  are  limited  to  the  E.  Kazakh  test 
site  (Case  6:  Table  2.2.3).  The  BEAM  parameters  produce  relatively  many  misclassified 
explosions,  with  one  exception  they  stem  from  the  E.  Kazakh  test  site. 

Among  the  earthquakes,  no.  26  in  lablc  2.2.2  appears  to  be  the  most  problematic,  since  it 
is  consistently  misclassified  both  by  BEAM  and  OGF  parameters,  except  for  one  maiginal 
rating  in  the  latter  case  (Case  6:  Table  2.2.3).  For  other  earthquakes  being  misclassified  or 
given  a  marginal  rating,  there  is  little  overlap  between  the  respective  OGF  and  BEAM  dis¬ 
criminators.  In  Fig.  2.2.8  number  of  events  arc  displayed  for  which  the  discriminators 
were  less  effective  or  failed.  We  would  specifically  comment  on  these  events  in  the  next 
.section. 

Discussion 

From  Table  2.2.3  we  have  that  the  poorest  performance  takes  place  for  the  combination  of 
all  events  and  all  discrimination  parameters  either  derived  from  OGF  or  BEAM  traces. 
Rather  surprisingly,  the  number  of  maiginal  events  seemingly  is  independent  of  event 
population  and  discrimination  parameters  used.  Regarding  the  mi.sclassified  NE  popula¬ 
tions,  most  of  the  events  here  are  non-Semipalatinsk  (E.  Kazakh),  the  exception  being 
event  22  (BEAM)  and  events  18  and  42  (OGF  -  Case  4:  Table  2.2.3).  Under  optimum  con¬ 
ditions  (Case  6;  Table  2.2.3),  event  42  is  labelled  marginal.  From  the  visual  inspections  of 
traces  (Fig.  9),  some  of  the  failures  are  rather  obvious,  but  this  is  not  the  case  for  NE  22 
(BEAM  -  Cases  1-3:  Table  2.2.3).  We  note  that  the  two  problematic  E.  Kazakh  events,  no 
22  and  42,  are  assigned  the  smallest  mb  values  (Table  2.2.1). 
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Also  some  of  the  presumed  earthquakes  are  problematic,  in  particular  EQ  11, 16  and  26 
(Cases  3  and  3:  Table  2.2.3).  The  BEAM  and  OGF  di.splays  in  Fig.2.2.8  give  the  visual 
impressions  that  EQ  1 1  and  26  exhibit  P  waves  rather  typical  for  explosions.  A  common 
feature  of  the  3  events  is  their  northernmost  latitude  of  47°  N,  while  corresponding  longi¬ 
tudes  are  around  83.3°,  89.7°  and  73.6°  E  (Table  2.2.2).  Also,  EQ  26  is  weak  (m5  4.3)  and 
besides  took  place  in  an  area  with  no  previous  seismic  activity.  The  ISC  bulletins  give  a 
normal  focal  depth  of  33  km.  and  that  Garm  is  the  only  USSR  reporting  station. 

Note,  most  discrimination  procedures  tend  to  favor  a  relatively  better  cltissification  perfor¬ 
mance  for  earthquakes  as  also  seen  in  Table  2.2.3.  The  rea.son  is  that  explosion  discrimi¬ 
nants  are  tied  to  relatively  weak  coda  excitation  and/or  P-signal  power  deficiency  in  low 
frctjucncy  bands.  I’or  decreasing  SNR  the  noise  contribution  would  be  more  relatively 
important  for  explosions  and  the  net  effect  is  that  all  weak  events  would  “look”  like  earth¬ 
quakes.  In  the  latter  case  it  seems  reasonable  to  select  for  learning  sets  weak  events  only 
and  then  take  advantage  of  the  OGF-technique  for  noise  suppression. 

From  the  above  we  have  that  the  performance  is  best  when  the  two  learning  populations 
are  drawn  from  roughly  the  same  area.  It  is  also  clear  that  the  di.scriminators  derived  from 
the  OGF-  traces  not  surprisingly  outpertbrm  the  coiresponding  BEAM  parameters.  This  is 
not  surpiising  for  the  simple  reason  that  the  signal-to-noise  ratio  is  significantly  better  on 
the  OGF  traces  (e.g.,  see  EQ  26  in  Fig.  2.2.8).  Not  easily  explainable  is  the  fact  that  the 
BEAM  parameters  are  weighted  in  favor  of  the  low  frequency  part  of  the  P  signal,  .slightly 
in  contrast  to  that  for  the  OGF  parameters. 

At  the  outset  of  this  study  the  “neural  network”  approach  was  not  contemplated,  and 
besides  Dy.sart  and  Pulli  (1990)  and  Dowlaet  al  (1990)  reported  that  the  Fisher  discrimi¬ 
nant  had  equivalent  performances.  The  reason  for  this  is  that  the  EQ  and  NE  populations 
are  e.ssentially  linearly  separated.  It  is  here  meant  that  the  feature  parameters  derived  from 
the  respective  event  populations  are  .separated  by  a  straight  line  in  a  two  parameter  case. 
We  would  here  add  that  our  approach  is  relatively  robust  for  small-sized  event  populations 
and  furthermore  that  the  probability  of  misclassification  introduced  here  is  helpful  in  judg¬ 
ing  discrimination  performances.  Further  details  are  given  in  Tsvang  et  al  (1992). 

Conclusion 

In  this  study  we  have  presented  a  comprehensive  seismic  source  discrimination  scheme 
for  NORESS-recorded  teleseismic  events,  comprising  44  nuclear  explosions  (NE)  and  35 
earthquakes  (EQ)  from  the  general  E.  Kazakh  area.  Major  results  were  as  follows: 

•  A  total  of  19  potential  discrimination  features  were  considered.  These  were  mainly 
tied  to  normalized  spectral  power  in  8  frequency  bands  for  both  P-  and  P-coda 
waves  of  durations  6  sec  and  30  sec,  respectively.  Also,  the  classical  P-complexity 
discriminant  was  incoiporated. 

•  The  most  effective  discrimination  features  were  the  complexity  one,  and  relatively 
low  P-  signal  frequencies.  The  best  classification  performance  was  obtained  by 
using  2  to  4  features  which  in  turn  were  selected  on  the  basis  of  the  extracted  mea¬ 
sure  of  maximum  misclassification  probability. 
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•  The  features  extracted  from  the  OGF  event  traces  had  a  slightly  better  performance 
than  those  derived  from  the  conventional  event  beam  traces.  The  corresponding 
optimum  feiture  subsets  were  somewhat  different. 

•  When  the  explosion  population  was  restricted  to  the  Semipalatinsk  area,  a  com¬ 
plete  event  classification  was  obtained  using  OGF-derived  features.  However,  3 
events  here  were  rated  of  marginal  significance. 

•  When  the  explosion  population  comprised  events  from  an  additional  7  testing 
areas,  the  classification  performances  decreased.  This  is  attributed  to  localized  site 
structures  and  upper  mantle  propagation  effects. 

•  A  remaining  problem  is  that  of  designing  discriminants  for  areas  where  available 
events  are  essentially  limited  to  either  earthquake  or  explosion  populations. 

S.L.  Tsvang,  V.I.  Pinsky,  MITPAN  Inst.,  Moscow 
E.S.  Husebye 


References 

Anderson,  T.  W.  (1958).  Introduction  to  multivariate  statistical  analysis.  Wiley,  New  York. 

Azen,  S.  P.,  L.  Breiman  and  W.  S.  Meisel  (1975).  Modem  approaches  to  data  analysis. 
Course  Notes.  Technology  Scivice  Corporation,  Santa  Monica,  CA. 

Oahlman,  0.  and  H.  Israelson  (1977).  Monitoring  of  Underground  Nuclear  Explosions. 
Elsevier,  New  York,  440  pp. 

Dainty,  A.  M.  ( 1990).  Studies  of  coda  using  array  and  three-component  processing, 
PAGEOPH,  132,221-244. 

Deev,  A.  D.  (1970).  Representation  of  statistics  of  discrete  analysis  and  asymptotic  expan¬ 
sion  for  .space  dimension  comparable  with  sample  size,  Dokl.  Akad.  Nauk,  v.l95,  N 
4,  USSR,  7.59-762. 

Dowla,  F.,  S.  Taylor  and  R.  Anderson  (1990).  Seismic  discrimination  with  artificial  neural 
networks:  preliminaiy  results  with  regional  spectral  data.  Bull.  Seism.  Soc.  Am.  80, 
1346-1373. 

Dysart,  P.  S.  and  J.  J,  Pulli  (1990).  Regional  seismic  event  classification  at  the  NORESS 
array:  .seismological  measurements  and  the  use  of  trained  neural  networks.  Bull. 
Seism.  Soc.  Am.  80.  1910-1933. 

Evemden,  J.  F.  ( 1 977).  Spectral  characteristics  of  the  P  codas  of  Eurasian  earthquakes  and 
e).plo.sion.s.  Bull.  Seism.  Soc.  Am.,  67, 1 153-1171, 


22 


Basic  Scistnolo^icai  Kcscaa'Ii 


Novcmilcr  1991 


Evemden,  J.  F.,  C.  B.  Archambeau  and  E.  Cranswick  (1986).  An  evaluation  of  .seismic 
decoupling  and  underground  nuclear  test  monitoring  using  high-frequency  data. 

Rev.  Geophys.,  24, 143-215. 

Ingate,  S.  F„  E.  S.  Husebye  and  A.  Christoffersson  (1985).  Regional  arrays  and  optimum 
data  processing  schemes.  Bull.  Seism.  Soc.  Am.,  75, 1155-1177, 

Kushnir,  A.  R,  E.  V.  Troitskiy  and  S.  L.  Tsvang  (1986).  Dependence  of  the  probability  of 
recognition  errors  of  the  size  of  the  sample  being  classified.  Computational  Seismol¬ 
ogy,  Alicnon  Press,  Inc.,  19, 77-82. 

Kushnir,  A.  R,  V.  M.  Lapshin,  V.  1.  Pinsky,  and  J.  Fyen  (1990).  Statistically  optimal  event 
detection  using  small  anay  data.  Bull.  Seism.  Soc.  Am.  80, 1934-1950, 

Raudis,  Sh.  Yu.  and  V.  S.  Pikyalis  ( 1 975).  Tabulation  of  the  expected  classification  error  of 
a  linear  discriminant  function  as  afunction  of  the  training  .sample.  In;  Statistical 
problems  of  control.  Institute  Fiziki  i  Matematiki  AN  LitSSR,  Vilnnius,  N  11, 81- 
120  (in  russian). 

Sandvin,  O.  and  D.  Tjostheim  ( 1978).  Multivariate  autoregressive  representation  of  seis¬ 
mic  P-wave  signals  with  application  to  short-period  discrimination.  Bull.  Seism. 

Soc.  Am.  68.  N3, 735-756. 

Tj0stheim,  D.  (1975).  Autoregressive  representation  of  seismic  P-wave  signals  with  an 
application  to  the  problem  of  short-period  discriminants.  Geophys.  J.  R.  astr.  Soc. 
43. 269-291. 

Tj0stheim.  D.  (1978).  Improved  seismic  discrimination  using  pattern  recognition,  Phys. 
Earth  Planet.  Int.,  16, 85-108. 

Tsvang,  S.  L..  A.  F.  Kushnir.  S.  P.  Starodubrovskaya,  N.  G.  Gamburtseva,  I.  A.  Ravich 
(1986).  Statistical  Analysis  of  Wave  Tails  in  Seismic  Sounding,  Izvestia,  Earth  Phys¬ 
ics.  22. 7. 540-552. 


Tsvang.  S.L.,  V.I.  Pinsky  and  E.S.  Husebye  ( 1992):  Enhanced  seismic  source  discrimina¬ 
tion  using  NORESS  recordings  from  Eurasian  events,  manuscript  submitted  for  pub¬ 
lication. 

Weber  K.,  A.  D.  Gvishiani,  P.  Gaudefrois,  A.  1.  Gorshkov,  A.  F.  Kushnir,  V.  F.  Pisarenko, 
A,  V.  V.  Trusov,  M.  L,  Tsvang  and  S.  L.  Tsvang  (1986).  Classification  of  high-seis¬ 
micity  Zones  in  the  Western  Alps,  Izvestia,  Earth  Physics,  22, 12, 965-976. 


23 


Basic  Scismologtcal  Research 


November  1991 


Presumed  Explosions 


No 

Date 

Region 

l.al. 

IvOn. 

Mb 

Ms 

1 

1985  02  10 

Ea.st.Kaz.  (5) 

49.87 

78.81 

5.9 

4.4 

2 

1984  12  28 

East.Kaz.  (5) 

49.86 

78.75 

6.0 

4.1 

3 

1988  0*1  14 

1-ast.Kaz.  (5) 

49.82 

78.79 

6.1 

4.6 

4 

1988  11  12 

Ea.st.Kaz.  (5) 

50.05 

78.99 

5.2 

— 

5 

1988  11  23 

East.Kaz.  (5) 

49.78 

78.14 

5.3 

— 

6 

1985  04  25 

l-asl-Kaz.  (5) 

49.92 

78.97 

5.9 

5.2 

7 

1987  02  26 

Ea.st.Kaz.  (5) 

49.81 

78.16 

5.4 

— 

8 

1987  03  12 

[^.st.Kaz.  (5) 

49.93 

78.78 

5.4 

3.9 

9 

1987  (M  03 

East.Kaz.  (5) 

49.90 

78.80 

6.2 

4.7 

10 

1987  04  17 

East.Kaz.  (5) 

49.85 

78.69 

6.0 

4.3 

11 

1987  04  19 

Ural.Mounl.  (8) 

60.78 

57.50 

4.5 

— 

12 

1987  05  06 

East.Kaz.  (5) 

49.77 

78.09 

5.5 

— 

13 

1987  06  05 

S.Sinkiang  (6) 

41.58 

88.75 

6.3 

4.4 

14 

1987  06  06 

East.Kaz.  (5) 

49.86 

78.14 

5.4 

— 

15 

1987  07  06 

Cent.Siberia  (3) 

61.49 

112.78 

5.2 

— 

16 

1987  07  17 

East.Kaz.  (5) 

49.78 

78.12 

5.8 

4.5 

17 

1987  07  24 

Cent.Siberia  (3) 

61.46 

112.72 

5.1 

— 

18 

1987  08  02 

East.Kaz.  (5) 

49.84 

78.88 

5.9 

3,8 

19 

1987  08  02 

N.Zcnilya  (1) 

73.31 

54.71 

5.8 

3.4 

20 

1987  08  12 

Cent.Siberia  (3) 

61.42 

112.71 

5.0 

... 

21 

1987  10  03 

Wcst.Kaz.  (4) 

47.63 

56.21 

5.2 

... 

22 

1987  10  16 

East.Kaz.  (5) 

49.78 

78.24 

4.6 

— 

23 

1987  11  15 

East.Kaz.  (5) 

49.87 

78.79 

6.0 

4.8 

24 

1987  12  13 

East.Kaz.  (5) 

49.95 

78.85 

6.1 

4.5 

25 

1987  12  20 

East.Kaz.  (5) 

49.75 

78.02 

4.8 

— 

26 

1987  12  27 

East.Kaz.  (5) 

49.83 

78.74 

6.1 

4.5 

27 

1988  02  06 

East.Kaz.  (5) 

49.86 

77.% 

4.8 

... 

28 

1988  02  13 

East.Kaz.  (5) 

49.92 

78.90 

6.0 

4.5 

29 

1988  04  03 

East.Kaz.  (5) 

49.88 

78.% 

5.9 

... 

30 

1988  05  04 

East.Kaz.  (5) 

49.91 

78.72 

6,2 

... 

31 

1988  05  07 

N.Zemlya  (1) 

73.35 

54.26 

5.5 

3.8 

32 

1988  08  22 

West.Siberia  (2) 

66.28 

78.55 

5.3 

... 

33 

1988  06  14 

East.Kaz.  (5) 

50.02 

78.98 

4.9 

4.1 

34 

1988  09  29 

S.Sinkiang  (6) 

41.82 

88.25 

4.8 

... 

35 

1988  10  18 

East.Kaz.  (5) 

49.86 

78.12 

4.9 

... 

36 

1988  12  04 

N.Zemlya  (1) 

73..36 

55.07 

5.9 

... 

37 

1988  12  17 

liast.Kaz.  (5) 

49.88 

78.92 

5.9 

4,.S 

38 

1989  01  22 

l-ast.Kaz.  (5) 

49.91 

78.86 

6.0 

4.5 

39 

1989  02  12 

lia.st.Kaz.  (5) 

49.89 

78.75 

5.9 

4.4 

40 

1989  02  17 

l^st.Kaz.  (5) 

49.87 

78.12 

5.0 

— 

41 

1985  06  30 

liast.Kaz.  (5) 

49.86 

78.69 

6.0 

4.2 

42 

1985  07  1 1 

East.Kaz.  (5) 

49.90 

78.80 

3.5 

... 

43 

1985  07  18 

Wcst.Russia  (7) 

65.97 

40.86 

5.0 

3.7 

44 

1985  07  25 

East.Kaz.  (5) 

49.89 

78.15 

5.0 

4.0 

Table  2.2.1.  Listing  of  the  presumed  underground  nuclear  explosions  recorded  at  NOR- 
ESS  and  used  in  our  analysis.  The  region  numbers  (in  brackets)  refer  to  the  specific  USSR 
test  sites  shown  in  Fig.  2.2.1.  The  focal  parameters  are  taken  from  the  PDE  and  ISC  bulle¬ 
tins,  Note  that  Lat.  and  Lon.  refer  to  latitude  (degrees  North)  and  longitude  (degrees  East), 
respectively.  Note,  Event  1 1  (Ural  Mountains)  is  the  second  explosion  on  19  April  with 
origin  time  at  04.04.55.6. 
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Presumed  luirlhquakcs 


No 

Date 

Ijit. 

Ijon. 

Depth 

Mb 

Ms 

1 

1984  11  02 

42.00 

84.06 

33 

4.5 

— 

2 

1985  04  16 

42.26 

82.24 

33(3) 

4.5 

3 

1986  05  30 

43.23 

87.82 

33 

4.6 

4 

1986  06  12 

43.67 

87.33 

33 

4.8 

5 

1986  07  03 

43.91 

84.69 

33 

4.4 

— 

6 

1986  07  17 

43.32 

77.85 

33 

4.5 

7 

198607  21 

44.66 

79.50 

33 

4.6 

... 

8 

198607  24 

43.79 

87.25 

33 

4.5 

9 

1986  10  04 

42.39 

84.63 

33 

4.0 

10 

1986  11  20 

42.04 

84.39 

50 

4.6 

... 

11 

1986  12  14 

47.31 

83.31 

33 

5.0 

... 

12 

1987  04  04 

42.45 

79.95 

33 

4.1 

4.1 

13 

1987  05  10 

44.28 

79.74 

33 

4.5 

14 

1987  05  26 

42.92 

78.06 

20 

4.6 

... 

15 

1987  08  22 

43.81 

85.29 

58 

4.4 

Ki 

1987  09  18 

47.01 

89.65 

33 

5.3 

4.8 

17 

1987  09  20 

42.91 

77.61 

41 

4.6 

4.2 

18 

1987  10  06 

43.43 

88.54 

32 

4.8 

4.0 

19 

1987  10  16 

44.20 

82.84 

56 

4.7 

... 

20 

1988  02  08 

43.73 

83.76 

10 

4.3 

... 

21 

1988  03  13 

42.18 

75.44 

33 

4.5 

... 

22 

198803  15 

42.21 

75.50 

33 

4.5 

... 

23 

1988  03  25 

44.70 

79.60 

33 

4.5 

... 

24 

198805  25 

42.01 

85.69 

22 

5.2 

25 

1988  06  17 

42.97 

77.50 

24 

5.3 

5.3 

26 

1988  09  27 

46.80 

73.59 

5 

4.3 

27 

1988  11  15 

42.01 

89.29 

33 

5.0 

4.3 

28 

1989  03  05 

42.51 

74.63 

33 

5.3 

4.1 

29 

1989  05  08 

44.83 

79.92 

33 

4.7 

... 

30 

1985  02  03 

42.06 

8434 

- 

4.5 

— 

31 

1985  03  24 

42.06 

77.62 

3 

43 

... 

32 

1985  06  02 

43i}0 

85.65 

21(45) 

4.9 

4.2 

33 

1985  07  16 

42.22 

82.36 

45(22) 

4.9 

4.2 

34 

1985  08  14 

42.13 

82.44 

33 

4.1 

... 

35 

1985  08  23 

42.65 

74.50 

3 

4.3 

... 

Table  2.2.2.  Listing  of  presumed  eanhquakes  recorded  at  NORESS  and  used  in  analysis. 
Epicenter  locations  are  shown  in  the  insert  in  Fig.  2.2.1.  Caption  otherwise  as  for  Table 
2.2.1. 
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DISCRIMINATION  RESULTS 


Case  Data  Populations  Class,  errors  Qass.  Misclass.  Marginal 

Rcclass.  J. knife  features  events  events 


1  BEAM  All  events  7 

2  BEAM  All  events 

3  BEAM  All  EQ 

Kazakh  NE 

4  OOF  All  events 

5  OGF  All  events 

6  OGF  All  EO 

Kazakh  NE 


9  All 

7  1,19,3,13 

1  3  19,13 

5  8  All 

3  5  19,10,15, 

13 

0  0  19,10 


NE:22,21,31,32, 

34,43; 

NE;11; 

EQ:7,12,32: 

EQ:  16, 17,29; 

NE:22,19,31,32, 

34,43; 

NE:36,42; 

EQ:17; 

EQ:7,25; 

NE:22; 

EQ:16,26; 

EQ:17; 

NE:18,19,32,34 

42,43; 

NE:21; 

EQ:11,26; 

EQ:7,16,17; 

NE:21,32,43 

EQ:11,26 

NE:19; 

" 

NE:25,30,42; 

EQ:7,26; 

Table  2.2.3.  Discrimination  results  for  the  Eurasian  events  listed  in  Tables  2.2.1  and  2.2.2, 
and  recorded  at  the  NORESS  array.  Classification  features  are  extracted  from  ordinary 
beam  (BEAM)  traces  and  optimunr  group  filtered  (OGF)  traces,  where  coherent  noise  has 
been  suppressed.  The  explosion  (NE)  population  was  divided  in  two  parts:  all  explosions 
or  only  those  at  the  E.  Kuzakh  test  site  (see  Fig.  2.2. 1 ).  Classification  (class.)  errors  are 
given  for  two  cases:  RecKassification  where  the  event  tested  was  pan  of  the  learning  popu¬ 
lation  and  Jack-knife,  where  test  event  and  learning  population  were  independent.  The 
“Reclass.”  results  are  always  somewhat  optimistic.  Classification  features  are  detailed  in 
Figs.  2.2.3  and  2.2.6.  Misclassified  events  and  Marginal  events  refer  to  the  numbering  in 
Table  2.2.1  and  2.2.2.  The  most  problematic  event  is  EQ  26,  that  is,  a  presumed  earth¬ 
quake  of  09/27/88  located  at  46.8N,  73.6E  (see  Fig.  2.2.1),  which  also  has  a  visual  appear¬ 
ance  of  being  an  explosion.  Waveforms  for  a  representative  number  of  events  being  either 
misclassified  or  rated  marginal  are  shown  in  Fig.  2.2.8. 
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Fig.  2.2.1.  Map  outlining  presumed  explosion  sites  and  presumed  earthquake  epicenters 
used  in  this  classification  study.  The  test  sites  are  N.  Zemlya  (1);  W.  Siberia  (2);  Cent. 
Siberia  (3);  W,  Kazakh  (4);  E.  Kazakh  (5);  S.  Sinkiang  (6);  W.  Russia  (7);  and  Ural  Moun¬ 
tains  (8).  Focal  parameter  details  in  Tables  2.2.1  and  2.2.2. 
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POWER  SPECTRA  (8  bands) 


►  ejim>  60.16 ;  vei*  2V79 
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0.01 


beam 


EQ:  15 
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0.0  10  20  XO  «0  SO  ■ 
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Fig.  2.2.2.  Example  of  the  relative  efficiency  of  the  Optimum  Group  Filtering  (OGF) 
scheme  of  Kushnir  et  al  (1990)  relative  to  conventional  array  beamforming.  The  OGF  fil¬ 
ter  suppresses  most  of  the  low  frequency  coherent  noise  as  clearly  seen  in  the  figure  to  the 
left.  The  event  analyzed  here  is  an  earthquake  -  no.  15  in  Table  2.2.2. 
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PARAMETERS 


P-phai*  (e  lie.) 


pifimtttfi  numbtf 

1 

2 

3 

4 

5 

6 

7 

8 

pirimittn 

0, 2-0.8 

0.8- 1.4 

1.4-2.0 

2.0-2.6 

2.6.3.2 

3.2.3.8 

3.8-4.4 

4.4-5.0 

POWER  SPECTRUM  (exomple) 


9 

10 

11 

12 

13 

14 

15 

16 

0.2.0.8 

CO 

o 

1.4-2.0 

2.0-2.6 

2.6-3.2 

3.2-3.8 

3.8-4.4 

4.4-5.0 

17  1 

18 

19 

FAmi.  (P-phase)l 

iFAmi,  (P-codaT 

Amax(P)/Amax(coda) 

Fig.  2.23.  The  19  feature  parameters  used  in  this  classification  study.  Center  frequencies 
of  the  spectral  bands  used  are  shown  to  the  left.  Note,  all  spectral  band  values  were  nor¬ 
malized  for  each  event  with  respect  to  Amax.  In  the  calculation,  the  logarithmic  feature 
values  were  used. 
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P-spectra 


explosions 
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t  2  3  4  5 
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Fig.  2.2.4.  Normalize  power  specira  (signal-noi.se)  for  all  eventa  used  in  analysis  (Table 
2.2.1  and  2.2.2).  “Peaks”  correspond  to  the  center  frequency  of  the  eight  spectral  bands 
shown  in  Fig.  2.2.2.  Seemingly,  the  largest  spectral  differences  between  the  explosion  and 
earthquake  populations  are  in  the  range  0.5  - 1.5  Hz. 
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Fig.  2.2.5.  Schematic  illustration  of  the  event  classification  procedure  used.  Further  details 
in  the  text. 
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paramctgrs  selection  (OGF) 


Ronked  porometers  (o’ol!  events;  4- -only  SemlpoL) 


PARAMEIERS  SELECTION  (BEAM) 


Ranked  porometefs  (o-otl  events;  4-only  Semlpat.) 

Fig.  2  J.6.  The  relative  petfomiance  of  feature  parameters  and  combinations  hereof  as  a 
function  of  the  Expected  Probability  of  Misclassification  (EPMC)  as  defined  in  the  text. 
CXjF  refers  to  optimum  group  filtered  traces  (removal  of  coherent  noise),  while  BEAM 
refers  to  ordinary  beam  traces.  Also,  the  upper  curves  refer  to  all  events,  while  the  lower 
ones  refer  to  the  subset  of  the  Semipalatinsk  or  E,  Kazakh  test  site  explosions  (Fig.  2.2.1). 
The  minimum  in  the  EPMC  parameter  implies  that  a  combination  of  2  or  4  classification 
features  is  optimum,  as  inclusion  of  additional  feature  parameters  actually  decreases  the 
performance  as  discussed  in  the  text. 
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BEAM  DISCRIMINATION 

35E0+44NE,all  porometcrs 


35E0+44NE,  optimal  porameters  (1,19,3,13} 


Only  E.Kazach  (32NE+35EQ),  optimol  paramelers  (19,13) 


Cv«nt$  numbars  ronktd  du#  <o  IDF  volu«s 


A-c.tfseUssHied  explosions  •'Oiisctessified  eofthquakot 

Fig.  2.2.7a.  Event  discrimination  with  feature  parameters  extracted  from  the  conventional 
beam  traces.  The  LDF  (linear  discrimination  function)  is  defined  in  the  text  and  Jack- 
Knife  refers  to  the  specific  way  of  testing  individual  events  independently  against  the 
learning  population.  The  events  are  ranked  relative  to  their  LDF  values,  and  thus  not  the 
way  done  in  Tables  2.2.1  and  2.2.2.  Events  with  LDF  values  within  +  0.2  units  (stippled 
lines)  are  somewhat  arbitrrjily  defined  as  marginal  events.  The  3  combinations  of  event 
populations  and  optimum  parameter  combinations  coincide  with  Case  1-3  in  Table  2.2.3. 
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35t04-44NE,oll  porometers 


35EQ+44NE,  opfimol  poramelers  (19,10,15,13) 


0  to  20  SO  40 
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Fig.  2.2.7b.  Event  discrimination  with  feature  parameters  extracted  from  the  OGF  traces 
(optimum  group  filtered).  Caption  otherwise  as  for  Fig.  2.2.7a. 
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Fig.  2.2.8.  Plot  of  OGF  and  conventional  earn  traces  for  raisclassified  and  marginal 
events  as  listed  in  Table  2.2.3.  The  focal  [  luameters  of  the  EQ  and  NE  events  are  shown. 
Most  noticeable  here  is  EQ26  which  exhibits  the  basic  feature  of  an  explosion  P  wave. 
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2.3  Ray>bascd  investigation  of  Lg  azimuth  anomalies  at  ARCESS 


Introduction 

In  a  previous  report  (Bostock  et  al,  1990)  we  have  discussed  the  application  of  the  ray 
tracing  technique  of  Bostock  and  Kennett  ( 1 990)  to  the  characterization  of  Lg  propagation 
to  the  NORESS  array.  We  have  considered  the  application  of  the  same  approach  of  repre¬ 
senting  the  major  features  of  the  Lg  phase  via  constructively  interfering  S  waves  multiply 
reflected  within  the  earth’s  aust  to  understanding  the  azimuth  anomalies  observed  at  the 
ARCESS  Jirray. 

Lg  azimuth  anomalies  observed  at  ARCESS 

Figure  2.3.1  summarizes  the  Lg  azimuth  anomalies  at  ARCESS  reported  during  the  first 
eight  weeks  of  operation  of  the  IMS  array  analysis  system  in  October  and  Novembef 
1989.  The  anomalies  are  also  given  in  Table  2.3.1.  The  azimuth  anomalies  are  plotted  cen¬ 
tered  at  the  source  point  scaled  by  relative  size  with  circles  representing  negative  anoma¬ 
lies  and  triangles  positive  azimuth  anomalies.  The  position  of  the  ARCESS  anay  is 
indicated  by  the  filled  triangle.  The  pattern  of  azimuth  anomalies  in  Fig.  2.3.1  shows  rapid 
changes  in  both  the  size  and  sign  of  the  anomalies  with  geographic  position.  There  is, 
however,  a  general  uend  for  the  anomaly  to  be  negative  for  events  along  the  Gulf  of  Esto¬ 
nia  and  more  positive  for  events  north  of  Leningrad. 

Scandinavian  crustal  model  and  ray  analysis 

Most  models  of  the  crustal  thickness  under  Fennoscandia  show  a  considerable  depression 
of  the  Moho  under  the  southern  part  of  Finland  and  this  would  be  expected  to  have  a  sig¬ 
nificant  influence  on  the  propagation  of  guided  waves  within  the  crust.  The  paths  from 
.sources  in  the  Gulf  of  Estonia  pass  close  to  this  feature  and  the  predicted  propagation  pat¬ 
terns  depend  quite  .strongly  on  the  detailed  shape  of  the  Moho. 

For  consistency  with  the  earlier  work  we  have  retained  the  .same  model  of  the  Moho  struc¬ 
ture  under  the  Fennoscandian  region  (Fig.  2.3.2)  and  combined  this  with  a  smoothed  ver¬ 
sion  of  the  surface  topography.  We  simplify  the  cnistal  model  to  a  single  layer  with 
constant  velocity  with  boundaries  defined  by  the  smoothed  topography  and  Moho  varia¬ 
tions.  For  each  source  positioit  we  consider  shooting  a  spray  of  rays,  with  a  fixed  phase 
velocity  at  the  source,  and  then  track  the  propagation  of  these  rays  through  the  structure.  A 
convenient  repre.sentation  of  the  propagation  process  is  provided  by  mapping  the  succes¬ 
sive  points  of  multiple  reflection  on  the  Moho  boundary.  For  a  waveguide  of  coastant 
thickne.ss  the.se  reflection  points  would  be  arcs  of  circles,  and  for  a  complex  three-dimen¬ 
sional  structure,  the  geometrical  patterns  of  reflection  points  provide  a  useful  measure  of 
the  distortion  of  the  wavefronts  spreading  from  the  source. 

We  have  found  that  shooting  rays  with  an  initial  phase  velocity  of  4.0  km/s  gives  a  good 
measure  of  the  propagation  process  for  Lg  in  our  highly  simplified  model.  It  would  be  po.s- 
sible  to  take  the  effects  of  sedimentary  structures  near  the  surface  into  account  by  more 
complex  ray  tracing,  but  little  additional  physical  insight  would  thereby  be  obtained.  In 
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Fig.  2.3.3,  we  display  the  ray  propagation  patterns  for  three  source  locations  representa¬ 
tive  of  the  clusters  of  events  represented  in  the  observed  anomaly  patterns  of  Fig.  2.3.1. 

For  a  source  in  the  Gulf  of  Estonia  (A  in  Fig.  2.3. 1 ),  Fig.  2.3.3a  displays  the  pattern  of 
multiple  reflections.  The  presence  of  the  depression  in  the  Moho  surface  has  the  effect  of 
diverting  energy  towards  the  west  and  also  of  modifying  the  paths  towards  the  ARCESS 
array  which  is  marked  by  a  triangle.  The  ray  pattern  varies  significantly  as  the  source  posi¬ 
tion  is  moved  along  the  Gulf  of  Estonia,  and  it  is  likely  that  the  large  variation  in  observed 
azimuths  from  this  region  is  influenced  by  the  interaction  of  the  Lg  phase  with  the  com¬ 
plex  Moho  structure  in  southern  Finland. 

For  a  source  to  the  north  of  Leningrad  (B  in  Fig.  2.3.1 ),  as  shown  in  Fig.  2.3.3b,  the  influ¬ 
ence  of  the  Moho  depression  is  more  subtle.  There  is  a  tendency  for  ray  paths  to  be  drawn 
towards  the  depression  because  of  the  influence  of  the  gradients  at  the  Moho.  The  result  is 
a  degree  of  def^ocusing  of  energy  travelling  in  the  general  direction  of  ARCESS.  Once 
again  the  propagation  patterns  depend  strongly  on  the  precise  location  of  the  source  rela¬ 
tive  to  the  Moho  depression  which  has  no  topographic  expression.  The  ray  diagrams  vary 
significantly  for  source  displacements  of  the  order  of  a  tenth  of  a  degree  which  again  is 
consistent  with  the  complexity  of  the  observations  from  this  area. 

A  more  northerly  source  (C  in  Fig.  2.3.1),  as  illustrated  in  Fig.  2.3.3c,  .shows  a  more  regu¬ 
lar  pattern  of  ray  propagation,  but  we  note  that  the  Moho  topography  in  northern  Finland 
has  led  to  a  disruption  of  the  reflected  wavefront  in  the  vicinity  of  the  ARCESS  array. 

Discussion 

The  azimuth  anomalies  calculated  from  this  simple  ray  tracing  scheme  are  given  in  Table 
2.3.2  and  arc  .seen  to  be  generally  consistent  with  the  trends  in  the  ob.served  anomalies  di.s- 
played  in  Fig.  2.3.1.  However,  the  ray  picture  for  Lg  is  probably  more  valuable  for  provid¬ 
ing  a  tool  to  investigate  the  .sensitivity  of  guided  wave  propagation  to  the  relative  position 
of  the  source  and  major  subsurface  features  such  as  the  significant  Moho  topography 
under  Finland. 

B.L.N.  Kennett,  M.G.  Bostock,  Research  School  of  Earth  Sciences, 

Australian  National  University 
S.  Mykkeitveit 
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Tabic  2.3.1.  Lg  azimuth  anomalies  at  ARCESS  from  IMS  analysis.  The  anomalies  are  rel¬ 
ative  to  “true"  azimuths  based  on  event  locations  in  the  monthly  bulletins  issued  by  the 
University  of  Helsinki. 


39 


c 


Spisniologicn)  Kc  m  .m  li 


iy,ovc;nib(er  i99l 

■  >.  ’n'  j. - 


Source 

Phase  Vel. 

Back  Azim. 

A  Azim 

A 

3.8 

175®.14 

-0.43® 

59.4  N 

4.C 

-1.40® 

27.1  E 

4.2 

-1.39® 

B 

3.8 

157®.81 

+3.75® 

61.6  N 

4.0 

+6,61® 

31.5  E 

4.2 

+7.99® 

C 

3.8 

151®.54 

-0,74® 

64.8  N 

4.0 

-0.73® 

30.5  E 
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+0.43® 

Tabl«  23.2.  Azimuth  anomalies  from  ray  tracing  for  sources  at  positions  A.  B  and  C  (Fig. 
2.3.1).  These  theoretical  anomalies  were  computed  for  phase  velocities  of  3.8, 4.0  and  4.2 
km/s.  The  sign  of  the  theoretical  anomalies  corresponds  to  the  same  convention  as  the 
observed  anomalies  in  Table  2.3. 1 . 
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Lg  Rayplot 
ARCESS  Analysis 
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Fig.  2.3.1.  The  figure  shows  Lg  azimuth  anomalies  at  ARCESS  reported  during  the  first 
eight  weeks  of  operation  of  IMS  in  October  and  November  1989.  Circles  represent  nega¬ 
tive  values  and  triangles  positive  values  of  the  anomalies.  The  position  of  ARCESS  is 
indicated  by  a  filled  triangle.  Also  marked  are  the  three  source  loca'  .ons  (A,B.C)  used  in 
the  ray  analysis  in  Fig.  2.3.3. 
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SCALES 

Horizontal  -  km  as  marked 
Vertical  -  km  below  msl 


Fig.  2.3.2.  Moho  relief  over  Scandinavia,  Finland  and  the  Baltic  states  after  Bannister  et  al 
(1991);  lighter  shades  indicate  increased  Moho  depths.  The  location  of  the  ARCESS  array 
is  marked  by  a  triangle. 
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Fig.  2.3.3a.  Ray  diagrams  for  a  source  located  at  59.4®N,  27.1®E  at  phase  velocities  3.8 
km/s  (top)  and  4.0  km/s  (bottom). 
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ARCESS  Analysis 
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ARCESS  Analysis 
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Fig.  2.3.3b.  Ray  diagrams  for  a  source  located  at  61.6®N,  31.5®E  at  phase  velocities  3.8 
km/s  (top)  and  4.0  km/s  (bottom). 
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ARCESS  Analysis 
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Fig.  2.3.3c.  Ray  diagrams  for  a  source  located  at  64,8®N,  30.5®E  at  phase  velocities  3.8 
km/s  (top)  and  4.0  km/s  (bottom). 
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2.4  Lg  as  a  yield  estimation  tool 

Introduction 

The  seismic  Lg  wave  propagates  in  the  continental  lithosphere  and  can  be  observed  from 
large  explosions  as  far  away  as  5000  km  in  shield  and  stable  platform  areas  (Nuttli,  1973; 
Baumgardt,  1985).  Lg  is  generally  considered  to  consist  of  a  superposition  of  many 
higher-mode  surface  waves  of  grup  velocities  near  3.5  km/s,  and  its  radiation  is  therefore 
expected  to  be  more  isotropic  than  that  of  P  waves.  Thus,  full  azimuthal  coverage  is  not 
essential  for  reliable  determination  of  Lg  magnitude.  Furthermore,  Lg  is  not  affected  by 
lateral  heterogeneities  in  the  upper  mantle,  which  can  produce  strong  focusing/defocusing 
effects  on  P-waves,  and  therefore  contribute  to  a  significant  uncertainty  in  P-based  mt, 
estimates. 

In  recent  years,  the  Lg  phase  has  emeiged  as  maybe  the  most  promising  tool  to  obtain  pre¬ 
cise  yield  e.stimates  by  seismic  means.  The  pioneering  work  by  Otto  Nuttli  of  Saint  Louis 
University  in  the  early  1970s  first  focused  attention  on  using  the  Lg  phase  fo»‘  seismic 
source  size  e.stimation,  and  Nuttli  developed  over  the  years  a  general  technique  for  mea¬ 
suring  Lg  magnitude  -  m|,(Lg)  -  along  any  source-receiver  path.  He  successfully  applied 
this  method  to  obtain  excellent  Lg-based  yield  estimates  for  NTS  explosions,  and  also 
showed  that  the  Lg  phase  could  be  used  for  teleseismic  yield  estimation  at  the  main  Soviet 
nuclear  test  sites. 

Nuttli  did  his  readings  exclusively  from  analog  .seismograms,  using  a  very  sophisticated 
interactive  analysis  method.  With  the  emergence  of  widespread  digital  recordings  in  the 
1980s,  the  focus  shifted  to  automatic  digital  processing,  using  in  particular  RMS  measure¬ 
ments  of  digitally  filtered  traces  (typically  0.6-3.0  Hz)  in  a  fixed  time  window  (typically  2 
minutes).  This  method  allowed  reliable  estimates  to made  even  at  very  low  SNR,  by 
using  a  noise  compensation  procedure.  Assuming  that  independent  yield  estimates  are 
available  for  calibration  purpo.ses,  this  method  appears  paiticularly  well  suited  for  TTBT 
monitoring  of  test  sites  of  limited  geographical  extent. 

Extensive  research  has  been  conducted  by  a  number  of  .scientists  over  the  past  decade 
investigating  various  aspects  of  Lg  generation,  propagation,  attenuation  and  modelling. 
While  recognizing  the  importance  and  impact  of  this  work,  it  will  lead  too  long  to  try  to 
cover  all  of  the.se  topics  in  a  brief  review  paper.  I  will  therefore  focus  on  the  observational 
aspects  of  Lg  as  a  tool  for  yield  estimation.  For  a  more  extensive  review  of  Lg  related 
developments,  reference  is  made  to  the  paper  by  Hansen,  Ringda!  end  Richards  (1990). 

In  this  paper  a  review  is  given  of  the  Lg  magnitude  estimation  methodology,  and  of  results 
obtained  with  regard  to  yield  e.stimation.  Results  are  briefly  summarized  for  each  of  three 
test  sites;  Nevada,  Semipalatinsk  and  Novaya  2!lemlya.  The  paper  is  concluded  with  a  short 
di.scussion  of  the  capabilities  and  limitation  of  Lg  as  a  yield  estimation  tool.  Some  very 
recent  results  are  included,  in  particular  regarding  the  u.se  of  Lg  for  Novaya  Zemlya  explo¬ 
sions. 

Lf>  results  at  NTS 

Much  pioneering  work  on  Lg  waves  was  done  in  the  1970s  and  1980s  by  Otto  Nuttli  of 
Saint  Louis  University.  Thus,  Nuttli  (1973)  proposed  that  “since  Lg  represents  a  higher- 
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mode  wave  traveling  with  minimum  group  velocity”  it  would  be  appropriate  to  relate 
amplitude  (A)  and  distance  (A)  via 

A  =  ((sinA)"’^]e"^^  (1) 

where  K  is  governed  by  the  source  strength,  and  y  is  the  coefficient  of  anelastic  attenua¬ 
tion.  The  quality  factor  Q  is  related  to  y  by  the  equation: 

y={nf)l(QU)  (2) 


where  U  is  the  group  velocity. 

With  the  goal  of  defining  a  magnitude  .scale,  ba.sed  on  Lg  observations  Nuttli  described  in 
detail  (Nuttli  1973, 1986a)  a  three-step  procedure  to  obtain  what  he  called  an  mb(Lg) 
value  for  an  earthquake  or  an  explosion  of  interest.  The  three  steps  were  as  follows: 

(i)  7  was  estimated  for  a  particular  source-receiver  path; 

(ii)  equation  ( 1 )  was  used  to  predict  an  amplitude  at  one  particular  distance  (he  chose  A 
corresponding  to  10  km  for  reference);  and 

(iii)  magnitude  was  assigned  via  the  formula 

mf,(Lg)  =  5.0  +  log[A(10km)/110] 
where  A(10  km)  is  the  amplitude,  in  microns,  resulting  from  (ii). 

For  22  nuclear  explosions  below  the  water  table  at  NTS,  Nuttli  (1986a)  showed  that  his 
mb(Lg)  values,  using  only  three  WWSSN  stations  in  the  western  U.S.,  were  remarkably 
well  correlated  with  the  logarithm  of  announced  yield.  He  proposed  a  best-fitting  line 
through  this  magnitude-yield  data,  from  which  magnitudes  had  a  standard  deviation  of 
only  about  O.OS.  Patton  (1988)  developed  computer-automated  measures  of  Lg  amplitude 
aiming  at  reporducing  Nuttli’s  NTS  results  (see  Fig.  2.4.1).  Patton  measured  Lg  ampli¬ 
tudes  from  digital  seismograms  in  two  ways  -  by  using  the  third-largest  peak  and  by  com¬ 
puting  the  RMS  amplitude  in  the  Lg  time  window  -  and  found  very  little  difference 
(around  0.01  magnitude  unit)  in  the  amount  of  scatter  about  regression  lines  using  the  two 
measures.  However,  he  found  that  standard  deviations  from  best-fitting  mb(Lg)  / 
log(yield)  relations  were  low,  0.07-0.08  magnitude  units,  only  if  explosions  were 
restricted  to  sub-regions  of  NTS  (Pahute  Mesa,  northern  Yucca  Flat,  southern  Yucca  Flat). 

Lg  results  at  Semipalatinsk 

Based  on  the  success  in  estimating  yields  for  NTS  explosions,  Nuttli  proceeded  to  apply 
the  same  magnitude-yield  relation,  together  with  Lg  signals  recorded  at  analog  WWSN 
stations  in  Eurasia,  to  estimate  the  yields  of  nuclear  explosions  at  the  main  Soviet  test  site 
(Shagan  River,  Semipalatinsk)  (Nuttli  1986b).  For  the  period  1978-1984,  after  the  150  kt 
Threshold  Test  Ban  Treaty  had  gone  into  effect,  his  yield  estimates  for  Shagan  River 
explosions  included  twenty  that  exceeded  the  treshold,  including  one  (1982  December  5) 
estimated  by  Nuttli  to  be  about  3(X)  kt.  While  acknowledging  the  pioneering  work 
involved  in  these  studies,  it  is  clear  that  the  gene  ally  low  signal-to-noise  ratios  and  the 
problematic  data  quality  of  these  analog  recordings  made  very  precise  measurements 
impoissible  to  attain,  a  fact  also  recognized  by  Nuttli  himself.  Also,  at  the  teleseismic  dis- 
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tances  for  which  Nuttli  had  Lg  data,  1900-4400  km.  yield  estimates  based  on  absolute 
measures  of  ground  motion  that  have  to  be  extrapolated  back  to  10  km  are  a  severe  test  of 
the  valicity  of  eq.  (1),  and,  even  if  eq.  (1)  is  appropriate,  are  very  sensitive  to  errors  in  y  by 
10-15%  would  result  in  yield  estimates  about  two  times  too  high. 

In  the  first  of  a  number  of  Lg  studies  undertaken  by  the  NORSAR  staff  during  the  1980s, 
Ringdal  (1983)  analyzed  digital  NORSAR  Lg  data  of  selected  Semipalatinsk  undeiground 
nuclear  explosions.  He  found  that  when  using  NORSAR  RMS  Lg  instead  of  P  waves 
recorded  at  NORSAR  to  estimate  source  size,  it  was  possible  effectively  to  eliminate  the 
magnitude  bias  relative  to  world-wide  mi,  observed  at  NORSAR  between  Degelen  and 
Shagan  River  explosions.  The  method  consisted  of  averaging  log  (RMS)  values  of  indi¬ 
vidual  NORSAR  channels,  filtered  in  a  band  0.6-3.0  Hz  in  order  to  enhance  Lg  signal-to- 
noise  ratio.  Ringdal  and  Hokland  (1987)  expanded  the  data  base,  and  introduce  a  noise 
compensation  procedure  to  improve  the  reliability  of  measurement  at  low  SNR  values. 
They  were  able  to  idenify  a  distinct  P-Lg  bias  between  the  northeast  and  southwest  por¬ 
tions  of  the  Shagan  River  test  site  (see  Fig.  2.4.2),  a  feature  that  was  confirmed  by  Ringdal 
and  Fyen  (1988)  using  Graefenberg  array  data.  Ringdal  and  Marshall  (1989)  combined  P 
and  Lg  based  source  size  estimators  to  estimate  the  yields  of  96  Shagan  River  explosions 
during  1%5-1988,  suing  data  on  the  cratering  explosion  15  January  1965  as  a  reference 
for  the  yield  calculations. 

Hansen,  Ringdal  and  Richards  (1990)  analyzed  available  data  from  stations  in  Qtina  and 
the  Soviet  Union,  and  found  that  RMS  Lg  of  Semipalatinsk  exploions  measured  at  these 
sUitions  showed  excellent  consistency  (see  Fig.  2.4.3).  They  concluded  that  for  explosions 
at  Semipalatinsk  with  good  signal-to-noise  ratio.  mi,(Lg)  may  be  estimated  at  single  sta¬ 
tions  with  an  accuracy  (one  standard  deviation)  of  about  0.03  magnitude  unit.  It  is  note¬ 
worthy  that  this  high  accuracy  was  consistently  obtained  for  a  variety  of  stations  at  veiy 
different  azimuths  and  distances,  even  though  the  basic  parameters  remained  exactly  as 
originally  proposed  for  NORSAR  recordings  (0.6-3.0  Hz  bandpass  filter,  RMS  window 
length  of  2  minutes,  centered  at  a  time  corresponding  to  a  group  velocity  of  3.5  km/s). 

A  possibility  to  compare  Lg  and  P  magnitudes  to  published  yields  for  Semipalatinsk 
explosions  has  recently  emeiged  with  the  recent  publication  by  Soviet  scientists  quoting 
yield  estimates  for  a  number  of  such  explosion  (Bocharov  et  al,  1989;  see  also  Vergino, 
1989a.b).  Since  these  explosions  were  all  conducted  prior  to  1973,  there  are  very  few 
available  high-quality  digital  records  of  the  events.  Nevertheless,  based  on  those  NOR¬ 
SAR  recordings  that  are  available,  the  correspondence  between  log  RMS  Lg  and  log  yield 
is  excellent  (Ringdal,  1990).  A  suite  of  analog  recordings  obtained  from  stations  within 
the  Soviet  Union  has  recently  been  digitized  and  made  available  as  part  of  bilateral  US- 
USSR  arrangements.  Israelson  (1991)  has  computed  RMS  Lg  of  these  recordings  and  has 
found  excellent  correspondence  with  the  yields  published  by  Bocharov  et  al  (1989). 

/./» results  at  Novayo  Zemlya 

No  yield  data  have  so  far  been  published  from  Soviet  sources  for  explosions  at  the  Novaya 
Zemlya  test  site.  Thus  an  assessment  of  the  Lg  potential  can  only  be  made  in  an  indirect 
way.  Ringdal  and  Fyen  (1991)  compared  NORSAR  and  Grafen^rg  RMS  Lg  for  Novaya 
Zemlya  explosions  using  the  s.mie  procedure  as  earlier  done  for  Semipalatinsk.  Fig.  2.4.4 
shows  the  propagation  paths  to  the  two  arrays  along  with  examples  of  recorded  Novaya 
ZemI)  a  explosions. 
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Fig.  2.4.5  .show.s  the  correspondence  between  RMS  Lg  at  NORSAR  and  Grafenbei)»,  and 
also  compares  various  other  source  size  estimators.  Note  that  world-wide  mj,  corresponds 
extremely  well  with  NORSAR  P-coda  magnitude,  but  correlates  very  poorly  with  either 
NORSAR  or  Grafenberg  RMS  Lg.  Note  also  that  the  RMS  Lg  correspondence  between 
Grafenberg  and  NORSAR  is  excellent,  with  an  orthogonal  standard  deviation  of  only 
0.035.  The  scatter  is  further  reduced  (to  0.025)  if  we  consider  only  events  with  at  least  5 
available  GRF  channels  (Ringdal  and  Fyen,  1991).  Thus,  we  obtain  the  same  clo.se  corre¬ 
spondence  between  Lg  observations  from  the.se  two  arrays  for  Novaya  Zemlya  explosions 
as  has  previously  been  observed  for  Semipalatinsk  events. 

With  the  current  lack  of  independently  obtained  calibration  data,  it  would  be  premature  to 
draw  any  firm  conclusions  as  to  the  relative  accuracy  of  RMS  Lg  in  estimating  yields  of 
these  explosions.  Nevertheless,  it  would  appear  that  the  close  grouping  in  RMS  Lg,  espe¬ 
cially  seen  for  the  NORSAR  data,  is  unlikely  to  be  a  coincidence.  It  would  seem  reason¬ 
able  to  conclude  that  this  group  of  explosions  has  very  nearly  the  same  yield,  in  spite  of 
the  divergence  in  mi,  estimates.  However,  additional  analysis,  in  particular  including 
available  Lg  data  from  Soviet  stations  for  this  event  set,  should  be  performed  in  order  to 
further  test  this  hypothesis.  Initial  results  from  processing  data  from  Soviet  stations  seem 
to  give  some  support  in  this  regard  (Israelson,  1991). 

Discussion 

Mo.st  of  the  evidence  for  the  “.stability”  of  Lg  magnitudes  is  indirect.  By  pairwise  com- 
pari.son  of  RMS  Lg  for  a  number  of  different  .source-receiver  paths,  it  has  been  demon¬ 
strated  that  single-station  RMS  Lg  can  be  u.sed  to  estimate  relative  magnitude  with  a 
remarkably  .small  scatter  (0.02-0.03  in  mj,  units  orthogonally).  This  was  first  shown  for 
Semipalatinsk  explosions,  but  has  recently  been  confirmed  also  for  explosions  at  Novaya 
Zemlya.  The  latter  observation  is  particularly  interesting,  since  many  of  the  paths  exhib¬ 
ited  significant  “Lg  blockage”  effects  (Baumgardt,  1990).  The  Lg  blockage  is  for  exam¬ 
ple  seen  on  NORSAR  recordings,  where  the  Lg  phase  is  relatively  weak  compared  to  P 
and  Sn. 

NORSAR  and  Grafenberg  tirray  measurements  of  Lg  waves  from  the  Shagan  River  area 
have  been  demonstrated  as  being  sufficiently  precise  to  allow  a  systematic  P-Lg  magni¬ 
tude  bias  to  be  identified  between  the  NE  and  SW  parts  of  that  site.  Available  yield  data 
from  Soviet  sources  indicate  that  Lg  magnitudes  show  belter  consistency  with  yield  than 
does  P-based  magnitudes  for  explosions  from  this  area.  The  reason  for  this  P-Lg  bias  has 
been  the  subject  of  much  discussion.  The  most  likely  explanation  appears  to  be  P-wave 
focusing/defocusing  effects  in  the  upper  mantle  underlying  the  source  region. 

A  heuristic  explanation  for  the  apparently  superior  stability  of  Lg  waves  compared  to  P 
waves  in  measuring  source  size  lies  in  the  difference  in  the  nature  of  the  sampling  of  the 
seismic  source  for  each  of  these  phases.  Teleseismic  P  waves  sample  only  a  very  limited 
ponion  of  the  focal  sphere,  and  are  susceptible  to  strong  focusing/defocusing  effects  in  the 
upper  mantle.  Lg  waves  are  composed  of  multiple  rays  that  sample  a  laiger  portion  of  the 
focal  sphere,  and  appear  to  have  less  significant  focusing/defocusing  effects  along  their 
propagation  paths.  In  a  sense,  Lg  waves  “let  the  Eanh  do  the  averaging”. 

Still,  there  are  many  uncertainties  remaining  with  regard  to  the  potential  of  Lg  as  a  yield 
estimator.  The  most  significant  question  mark  would  appear  to  be  the  effects  of  full  or  par- 
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tial  “Lg  blockage”.  As  noted  above,  such  blockage  does  not  seem  to  have  much  effect  on 
RMS  Lg  stability  if  the  blockage  structure  is  well  removed  from  the  source  region  (cf. 
NORSAR  Lg  blockage  for  Novaya  Zemlya  explosions).  This  is  reasonable  since  the 
blockage  effects  would  in  this  case  be  similar  for  all  the  events. 

However,  a  different  situation  arises  if  full  or  partial  “Lg  blockage”  structures  are  located 
within  or  in  the  immediate  neighborhood  of  the  test  site.  It  is  easy  to  see  that  the  .stability 
of  Lg  would  be  severely  affected  in  such  ca.ses.  In  this  connection,  it  is  noteworthy  that 
RMS  Lg  at  NTS  does  not  seem  to  show  the  same  “stability”  as  at  the  two  Soviet  test  sites 
(Richards,  personal  communication).  It  would  seem  likely  that  such  local  partial  “block* 
age”  effects  at  NTS  might  account  for  this  difference. 

Lg  waveforms  cannot  at  pre.sent  be  modelled  with  the  same  quality  of  fit  between  synthet¬ 
ics  and  data  that  has  been  attained  with  other  pha.ses,  and  many  aspects  of  Lg  generation 
and  propagation  characteristics  are  still  not  well  understood.  Among  topics  that  need  fur¬ 
ther  study  are  Lg  blockage  and  scattering  effects  caused  by  tectonic  heterogeneities  and 
the  effects  of  topography,  near  source  geology  and  depth  of  burial  on  Lg  excitation.  Much 
challenging  work  therefore  remains  to  be  done  in  this  field. 

F.  Ringdal 
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Fig.  2.4.2.  Map  of  the  Shagan  River  area  (top),  showing  the  surface  geology  (see  Leith, 
1987).  The  bottom  figure  illustrates  the  systematic  P-Lg  magnitude  bias  across  the  test 
site,  and  the  correction  with  the  observed  faults.  (After  Ringdal  and  Marshall,  1989.) 
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Fig.  2.4.3.  The  map  at  the  top  shows  the  location  of  stations  analyzed  in  the  study  by 
Hansen.  Ringdal  and  Richards  (1990)  (triangles)  relative  to  the  two  main  Soviet  test  sites. 
The  bottom  part,  figures  a)-d).  shows  the  excellent  RMS  Lg  correspondence  between  four 
of  the  stations  in  the  USSR/China  and  the  NORSAR  array  for  Shagan  River  explosion*. 
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RMS  Lg  Novaya  Zemlya 
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Fig.  2.4.5.  Plots  showing  the  correspondence  between  RMS  Lg  at  NORSAR  and  Grafen- 
berg  and  world  wide  nn,.  Note  that  neither  NORSAR  nor  GRF  Lg  correlates  well  with 
((a)  and  (b),  but  they  are  mutually  very  consistent  (d).  Also  note  the  excellent  correspon¬ 
dence  between  mb  and  NORSAR  P-coda  magnitude  (c). 
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2.5  Summaries  of  Quarterly  Technical  Reports  submitted 

During  FY91,  three  quarterly  technical  reports  were  submitted  on  this  contract.  The 
abstracts  of  these  papers  are  given  in  the  following. 

2.5.1  Scattering  of  regional  P„  by  Moho  topgraphy  -  T.  Kvaerna*  and  DJ. 
Doornbos^  NTNF/NORSAR,  Kjeller,  Norway,  and  ^  Institute  of 
Geophysics,  University  of  Oslo,  Norway) 

The  often  observed  relatively  large  amplitudes  in  the  later  part  of  the  P„  signal  cannot  be 
explained  with  the  traditional  interpretation  of  P„  in  1-D  crust-mantle  models.  To  deter¬ 
mine  the  cause  of  these  characteristics,  we  have  analyzed  in  some  detail  the  NORESS 
array  records  of  Pn  from  a  suite  of  quarry  blasts  in  S.W.  Norway.  Application  of  wide¬ 
band  frequency-wavenumber  analysis  to  these  records  confirms  that  there  is  a  slowness 
and  azimuth  anomaly  associated  with  the  dominant  part  of  the  wavetrain  and  that  it  is  con¬ 
fined  to  a  particular  frequency  range.  Moreover,  the  scattering  source  of  the  anomaly  is 
determined  to  be  at  Moho  depth,  which  is  consistent  with  the  concept  of  scattering  by 
topographic  relief.  We  demonstrate  the  viability  of  this  concept  by  means  of  numerical 
experiments,  showing  that  for  realistic  models  of  topography,  the  energy  flux  of  scattered 
P  can  dominate  the  specular  flux  (i.e.,  the  flux  in  the  direction  defined  by  the  ray  crossing 
a  plane  interface)  for  incidence  angles  approaching  the  critical  angle  of  P„.  Since  the 
effects  appear  to  be  systematic,  we  have  the  possibility  to  calibrate  the  Pn  parameters  for 
event  location  and  velocity  determination  purposes. 

2.5.2  Integrated  array  and  3-coniponent  processing  using  a  “microarray”  •• 

T.  Kvaerna  and  F.  Ringdal  (NTNF/NORSAR,  i^eller,  Norway) 

A  “miCToarray”  as  defined  in  this  paper  is  modeled  on  a  subgeometry  of  the  NORESS 
array  (Mykkeltveit  et  al,  1990)  and  comprises  a  3-component  center  seismometer  sur¬ 
rounded  by  3  closely  spaced  vertical-component  sensors  deployed  over  a  typical  aperture 
of  0.3  km.  Analysis  of  five  days  of  continuous  data  has  shown  that  such  a  system  com¬ 
bines  the  benefits  of  array  and  3-component  proce.ssing  in  providing  reliable  automatic 
detection,  phase  identification  and  location  of  weak  seismic  events  at  local  and  regional 
distances.  The  data  processing  has  comprised  a)  multiple-band  filtering,  b)  coherent  and 
incoherent  beamforming,  c)  STA/LTA  threshold  detection,  d)  broadband  frequenc>  wave- 
number  (f-k)  analysis  and  e)  automatic  pha.se  association  and  event  location.  Using  verti¬ 
cal  components  only,  broadband  f-k  array  analysis  enables  correct  phase  identification  (P- 
type  or  S-type  pha.se)  in  95  per  cent  of  the  cases,  and  gives  S-wave  azimuths  with  a  root- 
mean-square  (RMS)  deviation  of  13.9  degrees  from  the  estimates  of  the  full  NORESS 
array.  It  is  particularly  significant  that  the  small  array  eliminates  the  need  for  introducing 
particle  motion  models,  which  creates  ambiguities  in  3-component  jutalysis  of  secondary 
phases  when  interfering  SH  and  SV  phases  occur.  P-phase  azimuths  are  estimated  using 
integrated  airay  and  3-component  f-k  analysis,  and  have  an  RMS  deviation  relative  to 
NORESS  of  only  9.6  degrees.  Compared  to  the  full  NORESS  array,  the  P-wave  detection 
capability  is  good  for  events  with  epicenters  within  500  km  of  the  station,  but  for  greater 
distances  the  performance  is  significantly  reduced.  The  S-phase  detection  capability  is 
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enhanced  by  incoherent  beamforming  of  the  horizontal  channels,  and  approaches  that; of 
NORESS  at  all  distances.  A  considerable  reduction  in  the  detector  false  alapi  rate,  is 
achieved  by  imposing  constraints  on  the  estimates  of  apparent  velocity  obtained  from  the 
frk  analysis  before  accepting  a  detected  phase. 

Z.5.3  Diffraction  and  seismic  tomography  ••  DJ.  Doornbos  (Institute  of 
Geophysics,  University  of  Oslo,  Norway) 

Diffraction  tomography  is  formulated  in  such  a  way  that  the  data  (travel  time  —  ot  wave* 
form  perturbations)  are  related  to  the  medium  perturbations  through  the  sum  of  two  teus^, 
The  first  term  is  the  ray  integral  of  ordinaiy  tomography  and  involves  only  phase  perhvha- 
tions.  The  additional  diffraction  term  involves  both  phase  and  amplitude  peiturhadons.. 
The  diffraction  term  is  linear  in  the  gradients  of  the  velocity  perturbation  in  an  ^cpuittiiQ 
medium,  the  gradients  of  the  elastic  and  density  perturbations  in  an  elastic  mediua^t  and 
tfke  gradients  of  the  boundary  perturbations  the  wave  is  crossing.  This  formulatiprt  dte 
additional  advantage  that  unwanted  diffractions  from  the  nonphysical  boundajcy  of  the 
region  under  study  can  be  easily  removed.  Acoustic  scattering,  elastic  scattering,  and  scat* 
taring  by  boundary  perturbations  are  analyzed  separately.  Attention  is  paid  to  the  ade*' 
quacy  of  the  acoustic  approximation,  and  to  the  difference  between  perturbations  of  a 
boundary  level  (topography)  and  perturbations  of  boundary  conditions.  These  differences 
are  irrelevant  for  ordinary  seismic  tomography.  All  results  are  based  on  first-order  approx* 
imations  (Bom  or  Rytov),  as  is  the  case  for  other  published  methods  of  diffraction  tomog¬ 
raphy. 
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